Python scipy.interpolate.CubicSpline() Examples

def cubic_spline(orbit_data):
    '''
    Compute component wise cubic spline of points of input data

    Args:
        orbit_data (numpy array): array of orbit data points of the
        format [time, x, y, z]

    Returns:
        list: component wise cubic splines of orbit data points of the format [spline_x, spline_y, spline_z]
    '''
    time = orbit_data[:,:1]
    coordinates = list([orbit_data[:,1:2], orbit_data[:,2:3], orbit_data[:,3:4]])
    splines = list(map(lambda a:CubicSpline(time.ravel(),a.ravel()), coordinates))

    return splines 

def modbc(self):
        self.resorx = 0.0
        self.resory = 0.0

        for i in range(2, self.ni):
            xp = self.x(i, nj - 1)
            yp = self.y(i, nj - 1)
            xb = self.x(i, nj)
            yb = self.y(i, nj)
            ifail = 0
            call e02bcf(nicap7, xkn, cn, xb, 1, sn, ifail)
            interpolate.CubicSpline(x, y, axis=0, bc_type='not-a-knot', extrapolate=None)
            dydxb = min(sn(2), 1.d10)
            dydxb = max(sn(2), 1.d-10)
            if(sn(2).lt.0.0) dydxb = max(sn(2), -1.d10)
            if(sn(2).lt.0.0) dydxb = min(sn(2), -1.d-10)
            x(i, nj) = (yb-yp-xb*dydxb-xp/dydxb)/(-dydxb-1.0/dydxb)
            ifail = 0.0
            call e02bcf(nicap7, xkn, cn, x(i, nj), 1, sn, ifail)
            y(i, nj) = sn(1)
            resxb = abs(xb-x(i, nj))/xl
            resorx = resorx+resxb
            resyb = abs(yb-y(i, nj))/yl
            resory = resory+resyb 

def velocity_from_pseudo_enthalpy(self, pseudo_enthalpy, interp_type='CubicSpline'):
        """
        Returns the speed of sound in geometerized units in the
        neutron star at the specified pressure.

        Assumes the equation
        vs = c (de/dp)^{-1/2}

        :param pseudo_enthalpy (`float`): Dimensionless pseudo-enthalpy.
        :param interp_type (`str`): String specifying interpolation type.
                                    Current implementations are 'CubicSpline', 'linear'.

        :return v_s (`float`): Speed of sound at `pseudo-enthalpy` in geometerized units.
        """
        pressure = self.pressure_from_pseudo_enthalpy(pseudo_enthalpy, interp_type=interp_type)
        return self.dedp(pressure, interp_type=interp_type) ** -0.5 

def dedh(self, pseudo_enthalpy, rel_dh=1e-5, interp_type='CubicSpline'):
        """
        Value of [depsilon/dh](p)

        :param pseudo_enthalpy (`float`): Dimensionless pseudo-enthalpy.
        :param interp_type (`str`): String specifying interpolation type.
                                    Current implementations are 'CubicSpline', 'linear'.
        :param rel_dh (`float`): Relative step size in pseudo-enthalpy space.

        :return dedh (`float`): Derivative of energy-density with respect to pseudo-enthalpy
                                evaluated at `pseudo_enthalpy` in geometerized units.
        """

        # step size=fraction of value
        dh = pseudo_enthalpy * rel_dh

        eps_upper = self.energy_density_from_pseudo_enthalpy(pseudo_enthalpy + dh, interp_type=interp_type)
        eps_lower = self.energy_density_from_pseudo_enthalpy(pseudo_enthalpy - dh, interp_type=interp_type)

        return (eps_upper - eps_lower) / (2. * dh) 

def __init__(self, p_1, p_2, duration=1.0):
        self.cs = CubicSpline(np.array([0.0, duration]), np.array([p_1, p_2]), bc_type='clamped') 

def k2_from_mass(self, m):
        """
        :param m: mass of neutron star in solar masses.
        :return: dimensionless second tidal love number.
        """
        f = CubicSpline(self.mass, self.k2love_number, bc_type='natural', extrapolate=True)

        m_geom = m * MSUN_SI * G_SI / C_SI ** 2.
        return f(m_geom) 

def radius_from_mass(self, m):
        """
        :param m: mass of neutron star in solar masses
        :return: radius of neutron star in meters
        """
        f = CubicSpline(self.mass, self.radius, bc_type='natural', extrapolate=True)

        mass_converted_to_geom = m * MSUN_SI * G_SI / C_SI ** 2.
        return f(mass_converted_to_geom) 

def pseudo_enthalpy_from_energy_density(self, energy_density, interp_type='CubicSpline'):
        """
        Find h(epsilon)
        as in lalsimulation, return h = K * e**(2./3.) below min enthalpy

        :param energy_density (`float`): energy-density in geometerized units.
        :param interp_type (`str`): String specifying interpolation type.
                                    Current implementations are 'CubicSpline', 'linear'.

        :return pseudo_enthalpy (`float`): Dimensionless pseudo-enthalpy.
        """
        energy_density = np.atleast_1d(energy_density)
        pseudo_enthalpy_returned = np.zeros(energy_density.size)
        indices_less_than_min = np.nonzero(energy_density < self.minimum_energy_density)
        indices_greater_than_min = np.nonzero(energy_density >= self.minimum_energy_density)

        pseudo_enthalpy_returned[indices_less_than_min] = 10 ** (np.log10(self.pseudo_enthalpy[0]) + (2. / 3.) *
                                                                 (np.log10(energy_density[indices_less_than_min]) -
                                                                  np.log10(self.energy_density[0])))

        if interp_type == 'CubicSpline':
            x = np.log10(energy_density[indices_greater_than_min])
            pseudo_enthalpy_returned[indices_greater_than_min] = 10**self.interp_pseudo_enthalpy_from_energy_density(x)
        elif interp_type == 'linear':
            pseudo_enthalpy_returned[indices_greater_than_min] = np.interp(energy_density[indices_greater_than_min],
                                                                           self.energy_density,
                                                                           self.pseudo_enthalpy)
        else:
            raise ValueError('Interpolation scheme must be linear or CubicSpline')

        if pseudo_enthalpy_returned.size == 1:
            return pseudo_enthalpy_returned[0]
        else:
            return pseudo_enthalpy_returned 

def energy_density_from_pseudo_enthalpy(self, pseudo_enthalpy, interp_type='CubicSpline'):
        """
        Find energy_density_from_pseudo_enthalpy(pseudo_enthalpy)
        as in lalsimulation, return e = K * h**(3./2.) below min enthalpy

        :param pseudo_enthalpy (`float`): Dimensionless pseudo-enthalpy.
        :param interp_type (`str`): String specifying interpolation type.
                                    Current implementations are 'CubicSpline', 'linear'.

        :return energy_density (`float`): energy-density in geometerized units.
        """
        pseudo_enthalpy = np.atleast_1d(pseudo_enthalpy)
        energy_returned = np.zeros(pseudo_enthalpy.size)
        indices_less_than_min = np.nonzero(pseudo_enthalpy < self.minimum_pseudo_enthalpy)
        indices_greater_than_min = np.nonzero(pseudo_enthalpy >= self.minimum_pseudo_enthalpy)

        energy_returned[indices_less_than_min] = 10 ** (np.log10(self.energy_density[0]) +
                                                        1.5 * (np.log10(pseudo_enthalpy[indices_less_than_min]) -
                                                               np.log10(self.pseudo_enthalpy[0])))
        if interp_type == 'CubicSpline':
            x = np.log10(pseudo_enthalpy[indices_greater_than_min])
            energy_returned[indices_greater_than_min] = 10 ** self.interp_energy_density_from_pseudo_enthalpy(x)
        elif interp_type == 'linear':
            energy_returned[indices_greater_than_min] = np.interp(pseudo_enthalpy[indices_greater_than_min],
                                                                  self.pseudo_enthalpy,
                                                                  self.energy_density)
        else:
            raise ValueError('Interpolation scheme must be linear or CubicSpline')

        if energy_returned.size == 1:
            return energy_returned[0]
        else:
            return energy_returned 

def pressure_from_pseudo_enthalpy(self, pseudo_enthalpy, interp_type='CubicSpline'):
        """
        Find p(h)
        as in lalsimulation, return p = K * h**(5./2.) below min enthalpy

        :param pseudo_enthalpy (`float`): Dimensionless pseudo-enthalpy.
        :interp_type (`str`): String specifying interpolation type.
                              Current implementations are 'CubicSpline', 'linear'.

        :return pressure (`float`): pressure in geometerized units.
        """
        pseudo_enthalpy = np.atleast_1d(pseudo_enthalpy)
        pressure_returned = np.zeros(pseudo_enthalpy.size)
        indices_less_than_min = np.nonzero(pseudo_enthalpy < self.minimum_pseudo_enthalpy)
        indices_greater_than_min = np.nonzero(pseudo_enthalpy >= self.minimum_pseudo_enthalpy)

        pressure_returned[indices_less_than_min] = 10. ** (np.log10(self.pressure[0]) +
                                                           2.5 * (np.log10(pseudo_enthalpy[indices_less_than_min]) -
                                                                  np.log10(self.pseudo_enthalpy[0])))

        if interp_type == 'CubicSpline':
            pressure_returned[indices_greater_than_min] = (
                10. ** self.interp_pressure_from_pseudo_enthalpy(np.log10(pseudo_enthalpy[indices_greater_than_min])))
        elif interp_type == 'linear':
            pressure_returned[indices_greater_than_min] = np.interp(pseudo_enthalpy[indices_greater_than_min],
                                                                    self.pseudo_enthalpy,
                                                                    self.pressure)
        else:
            raise ValueError('Interpolation scheme must be linear or CubicSpline')

        if pressure_returned.size == 1:
            return pressure_returned[0]
        else:
            return pressure_returned 

def energy_from_pressure(self, pressure, interp_type='CubicSpline'):
        """
        Find value of energy_from_pressure
        as in lalsimulation, return e = K * p**(3./5.) below min pressure

        :param pressure (`float`): pressure in geometerized units.
        :param interp_type (`str`): String specifying which interpolation type to use.
                                    Currently implemented: 'CubicSpline', 'linear'.

        :param energy_density (`float`): energy-density in geometerized units.
        """
        pressure = np.atleast_1d(pressure)
        energy_returned = np.zeros(pressure.size)
        indices_less_than_min = np.nonzero(pressure < self.minimum_pressure)
        indices_greater_than_min = np.nonzero(pressure >= self.minimum_pressure)

        # We do this special for less than min pressure
        energy_returned[indices_less_than_min] = 10 ** (np.log10(self.energy_density[0]) +
                                                        (3. / 5.) * (np.log10(pressure[indices_less_than_min]) -
                                                                     np.log10(self.pressure[0])))

        if interp_type == 'CubicSpline':
            energy_returned[indices_greater_than_min] = (
                10. ** self.interp_energy_density_from_pressure(np.log10(pressure[indices_greater_than_min])))
        elif interp_type == 'linear':
            energy_returned[indices_greater_than_min] = np.interp(pressure[indices_greater_than_min],
                                                                  self.pressure,
                                                                  self.energy_density)
        else:
            raise ValueError('Interpolation scheme must be linear or CubicSpline')

        if energy_returned.size == 1:
            return energy_returned[0]
        else:
            return energy_returned 

def interpolation(data, ori_t_unit=0.02, tar_t_unit=0.01):
    assert(len(data.shape)==2)

    ori_x = np.arange(len(data))
    tar_x = np.arange(0, len(data), tar_t_unit/ori_t_unit)
    func = CubicSpline(ori_x, data, axis=0)
    return func(tar_x) 

def __init__(self, gate_poses):
        self.gates = gate_poses

        self.n_gates = np.size(gate_poses, 0)
        positions = np.array([pose.position.to_numpy_array() for pose in gate_poses])
        dists = np.linalg.norm(positions[1:, :] - positions[:-1, :], axis=1)
        self.arc_length = np.zeros(shape=self.n_gates)
        self.arc_length[1:] = np.cumsum(dists)

        # tangents from quaternion
        # by rotating default gate direction with quaternion
        self.tangents = np.zeros(shape=(self.n_gates, 3))
        for i, pose in enumerate(gate_poses):
            self.tangents[i, :] = rotate_vector(pose.orientation, gate_facing_vector).to_numpy_array()
        self.track_spline = CubicHermiteSpline(self.arc_length, positions, self.tangents, axis=0)

        # gate width to track (half) width
        gate_widths = [gate_dimensions[0] / 2.0 for gate in gate_poses]
        gate_heights = [gate_dimensions[1] / 2.0 for gate in gate_poses]

        self.track_width_spline = CubicSpline(self.arc_length, gate_widths, axis=0)
        self.track_height_spline = CubicSpline(self.arc_length, gate_heights, axis=0)

        # sample 2048 points, the 2048 are arbitrary and should really be a parameter
        taus = np.linspace(self.arc_length[0], self.arc_length[-1], 2**12)

        self.track_centers = self.track_spline(taus)
        self.track_tangents = self.track_spline.derivative(nu=1)(taus)
        self.track_tangents /= np.linalg.norm(self.track_tangents, axis=1)[:, np.newaxis]
        self.track_normals = np.zeros_like(self.track_tangents)
        self.track_normals[:, 0] = -self.track_tangents[:, 1]
        self.track_normals[:, 1] = self.track_tangents[:, 0]
        self.track_normals /= np.linalg.norm(self.track_normals, axis=1)[:, np.newaxis]

        self.track_widths = self.track_width_spline(taus)
        self.track_heights = self.track_height_spline(taus) 

def test_cubic_bc_natural():
    np.random.seed(1234)
    x = np.linspace(0, 5, 20)
    xi = np.linspace(0, 5, 100)
    y = np.sin(x) + np.random.randn(x.size) * 0.3

    cs = CubicSpline(x, y, bc_type='natural')
    ss = csaps.CubicSmoothingSpline(x, y, smooth=1.0)

    y_cs = cs(xi)
    y_ss = ss(xi)

    assert cs.c == pytest.approx(ss.spline.c)
    assert y_cs == pytest.approx(y_ss) 

def __init__(self, points, duration=1., bc_type='clamped'):
        n_points = points.shape[0]
        self._duration = duration
        self._cs = CubicSpline(np.linspace(0, duration, n_points), points, bc_type=bc_type) 

def stationary_rotation(dt, lat, alt, Cnb, Cbs=None):
    """Simulate readings on a stationary bench.

    Parameters
    ----------
    dt : float
        Time step.
    lat : float
        Latitude of the place.
    alt : float
        Altitude of the place.
    Cnb : ndarray, shape (n_points, 3, 3)
        Body attitude matrix.
    Cbs : ndarray with shape (3, 3) or (n_points, 3, 3) or None
        Sensor assembly attitude matrix relative to the body axes. If None,
        (default) identity attitude is assumed.

    Returns
    -------
    gyro, accel : ndarray, shape (n_points - 1, 3)
        Gyro and accelerometer readings.
    """
    n_points = Cnb.shape[0]

    time = dt * np.arange(n_points)
    lon_inertial = np.rad2deg(earth.RATE) * time
    lat = np.full_like(lon_inertial, lat)
    Cin = dcm.from_llw(lat, lon_inertial)

    R = transform.lla_to_ecef(lat, lon_inertial, alt)
    v_s = CubicSpline(time, R).derivative()

    if Cbs is None:
        Cns = Cnb
    else:
        Cns = util.mm_prod(Cnb, Cbs)

    Cis = util.mm_prod(Cin, Cns)
    Cib_spline = RotationSpline(time, Rotation.from_matrix(Cis))
    a = Cib_spline.interpolator.c[2]
    b = Cib_spline.interpolator.c[1]
    c = Cib_spline.interpolator.c[0]

    g = earth.gravitation_ecef(lat, lon_inertial, alt)
    a_s = v_s.derivative()
    d = a_s.c[1] - g[:-1]
    e = a_s.c[0] - np.diff(g, axis=0) / dt

    d = util.mv_prod(Cis[:-1], d, at=True)
    e = util.mv_prod(Cis[:-1], e, at=True)

    gyros, accels = _compute_readings(dt, a, b, c, d, e)

    return gyros, accels 

def _augment_core(
        self, X: np.ndarray, Y: Optional[np.ndarray]
    ) -> Tuple[np.ndarray, Optional[np.ndarray]]:
        N, T, C = X.shape
        rand = np.random.RandomState(self.seed)

        if isinstance(self.n_drift_points, int):
            n_drift_points = set([self.n_drift_points])
        else:
            n_drift_points = set(self.n_drift_points)

        ind = rand.choice(
            len(n_drift_points), N * (C if self.per_channel else 1)
        )  # map series to n_drift_points

        drift = np.zeros((N * (C if self.per_channel else 1), T))
        for i, n in enumerate(n_drift_points):
            if not (ind == i).any():
                continue
            anchors = np.cumsum(
                rand.normal(size=((ind == i).sum(), n + 2)), axis=1
            )  # type: np.ndarray
            interpFuncs = CubicSpline(
                np.linspace(0, T, n + 2), anchors, axis=1
            )  # type: Callable
            drift[ind == i, :] = interpFuncs(np.arange(T))
        drift = drift.reshape((N, -1, T)).swapaxes(1, 2)
        drift = drift - drift[:, 0, :].reshape(N, 1, -1)
        drift = drift / abs(drift).max(axis=1, keepdims=True)
        if isinstance(self.max_drift, (float, int)):
            drift = drift * self.max_drift
        else:
            drift = drift * rand.uniform(
                low=self.max_drift[0],
                high=self.max_drift[1],
                size=(N, 1, C if self.per_channel else 1),
            )

        if self.kind == "additive":
            if self.normalize:
                X_aug = X + drift * (
                    X.max(axis=1, keepdims=True) - X.min(axis=1, keepdims=True)
                )
            else:
                X_aug = X + drift
        else:
            X_aug = X * (1 + drift)

        if Y is not None:
            Y_aug = Y.copy()
        else:
            Y_aug = None

        return X_aug, Y_aug 

def _interpolate_packed_pixels(dataset, max_nans_interpolation):
        given_pos = np.arange(5, 409, 8)
        new_pos = np.arange(1, 410)

        lat_in = np.deg2rad(dataset["lat"].values)
        lon_in = np.deg2rad(dataset["lon"].values)

        # We cannot define given positions for each scanline, but we have to
        # set them for all equally. Hence, we skip every scan position of all
        # scan lines even if only one contains a NaN value:
        nan_scnpos = \
            np.isnan(lat_in).sum(axis=0) + np.isnan(lon_in).sum(axis=0)
        valid_pos = nan_scnpos == 0

        if valid_pos.sum() < 52 - max_nans_interpolation:
            raise ValueError(
                "Too many NaNs in latitude and longitude of this AVHRR file. "
                "Cannot guarantee a good interpolation!"
            )

        # Filter NaNs because CubicSpline cannot handle it:
        lat_in = lat_in[:, valid_pos]
        lon_in = lon_in[:, valid_pos]
        given_pos = given_pos[valid_pos]

        x_in = np.cos(lon_in) * np.cos(lat_in)
        y_in = np.sin(lon_in) * np.cos(lat_in)
        z_in = np.sin(lat_in)

        xf = CubicSpline(given_pos, x_in, axis=1, extrapolate=True)(new_pos)
        yf = CubicSpline(given_pos, y_in, axis=1, extrapolate=True)(new_pos)
        zf = CubicSpline(given_pos, z_in, axis=1, extrapolate=True)(new_pos)
        lon = np.rad2deg(np.arctan2(yf, xf))
        lat = np.rad2deg(np.arctan2(zf, np.sqrt(xf ** 2 + yf ** 2)))

        dataset["lat"] = ("scnline", "scnpos"), lat
        dataset["lon"] = ("scnline", "scnpos"), lon

        # The other packed variables will be simply padded:
        for var_name, var in dataset.data_vars.items():
            if "packed_pixels" not in var.dims:
                continue

            nan_scnpos = np.isnan(var).sum(axis=0)
            valid_pos = nan_scnpos == 0
            given_pos = np.arange(5, 409, 8)[valid_pos]

            dataset[var_name] = xr.DataArray(
                CubicSpline(
                    given_pos, var.values[:, valid_pos], axis=1,
                    extrapolate=True)(new_pos),
                dims=("scnline", "scnpos")
            )