Python scipy.interpolate.UnivariateSpline() Examples

def compute_depth_gravity_profiles(pressures, densities, surface_gravity, outer_radius):
    gravity = [surface_gravity] * len(pressures) # starting guess
    n_gravity_iterations = 5
    for i in range(n_gravity_iterations):    
        # Integrate the hydrostatic equation
        # Make a spline fit of densities as a function of pressures
        rhofunc = UnivariateSpline(pressures, densities)
        # Make a spline fit of gravity as a function of depth
        gfunc = UnivariateSpline(pressures, gravity)
        
        # integrate the hydrostatic equation
        depths = np.ravel(odeint((lambda p, x: 1./(gfunc(x) * rhofunc(x))), 0.0, pressures))
        
        radii = outer_radius - depths
        
        rhofunc = UnivariateSpline(radii[::-1], densities[::-1])
        poisson = lambda p, x: 4.0 * np.pi * burnman.constants.G * rhofunc(x) * x * x
        gravity = np.ravel(odeint(poisson, surface_gravity*radii[0]*radii[0], radii))
        gravity = gravity / radii / radii
    return depths, gravity


# BEGIN USER INPUTS

# Declare the rock we want to use 

def set_table(self, n=100, dx=None):
        r'''Method to set an interpolation table of liquids levels versus
        volumes in the tank, for a fully defined tank. Normally run by the
        h_from_V method, this may be run prior to its use with a custom
        specification. Either the number of points on the table, or the
        vertical distance between steps may be specified.

        Parameters
        ----------
        n : float, optional
            Number of points in the interpolation table, [-]
        dx : float, optional
            Vertical distance between steps in the interpolation table, [m]
        '''
        if dx:
            self.heights = linspace(0.0, self.h_max, int(self.h_max/dx)+1)
        else:
            self.heights = linspace(0.0, self.h_max, n)
        self.volumes = [self.V_from_h(h) for h in self.heights]
        from scipy.interpolate import UnivariateSpline
        self.interp_h_from_V = UnivariateSpline(self.volumes, self.heights, ext=3, s=0.0)
        self.table = True 

def calculate_service_battery_currents_v1(
        service_battery_capacity, times, service_battery_state_of_charges):
    """
    Calculate the service battery current vector [A].

    :param service_battery_capacity:
        Service battery capacity [Ah].
    :type service_battery_capacity: float

    :param times:
        Time vector [s].
    :type times: numpy.array

    :param service_battery_state_of_charges:
        State of charge of the service battery [%].
    :type service_battery_state_of_charges: numpy.array

    :return:
        Service battery current vector [A].
    :rtype: numpy.array
    """
    from scipy.interpolate import UnivariateSpline as Spline
    soc = service_battery_state_of_charges
    ib = Spline(times, soc, w=np.tile(10, times.shape[0])).derivative()(times)
    return ib * (service_battery_capacity * 36.0) 

def fit_policy_function(self,PF):
        '''
        Fits the policy functions PF using the points xgrid using UnivariateSpline
        '''
        xgrid,S = self.xgrid,self.S

        Vf,cf,nf,xprimef = {},{},{},{}
        for s in range(S):
            PFvec = np.vstack(map(lambda x:PF(x,s),xgrid))
            Vf[s] = UnivariateSpline(xgrid,PFvec[:,0],s=0)
            cf[s] = UnivariateSpline(xgrid,PFvec[:,1],s=0,k=1)
            nf[s] = UnivariateSpline(xgrid,PFvec[:,2],s=0,k=1)
            for sprime in range(S):
                xprimef[s,sprime] = UnivariateSpline(xgrid,PFvec[:,3+sprime],s=0,k=1)

        return Vf,[cf,nf,xprimef] 

def _compute_gravity(self, density, gravity_bottom):
        """
        Computes the gravity of a layer
        Used by _evaluate_eos()
        """
        # Create a spline fit of density as a function of radius
        rhofunc = UnivariateSpline(self.radii, density)
        # Numerically integrate Poisson's equation

        def poisson(p, x): return 4.0 * np.pi * \
            constants.G * rhofunc(x) * x * x
        grav = np.ravel(
            odeint( poisson, gravity_bottom *
                self.radii[0] * self.radii[0], self.radii))
        
        if self.radii[0] == 0:
            grav[0] = 0
            grav[1:] = grav[1:] / self.radii[1:] / self.radii[1:]
        else:
            grav[:] = grav[:] / self.radii[:] / self.radii[:]
        return grav 

def fit_policy_function(self, PF):
        '''
        Fits the policy functions PF using the points xgrid using UnivariateSpline
        '''
        xgrid, S = self.xgrid, self.S

        Vf, cf, nf, xprimef = {}, {}, {}, {}
        for s in range(S):
            PFvec = np.vstack(map(lambda x: PF(x, s), xgrid))
            Vf[s] = UnivariateSpline(xgrid, PFvec[:, 0], s=0)
            cf[s] = UnivariateSpline(xgrid, PFvec[:, 1], s=0, k=1)
            nf[s] = UnivariateSpline(xgrid, PFvec[:, 2], s=0, k=1)
            for sprime in range(S):
                xprimef[s, sprime] = UnivariateSpline(
                    xgrid, PFvec[:, 3 + sprime], s=0, k=1)

        return Vf, [cf, nf, xprimef] 

def linewidth(f, a):
    """Calculate the full-width at half maximum (FWHM) of an absorption line.

    Parameters:
        f (ndarray): Frequency grid.
        a (ndarray): Line properties
            (e.g. absorption coefficients or cross-sections).

    Returns:
        float: Linewidth.

    Examples:
        >>> f = np.linspace(0, np.pi, 100)
        >>> a = np.sin(f)**2
        >>> linewidth(f, a)
        1.571048056449009
    """
    s = interpolate.UnivariateSpline(f, a - np.max(a)/2, s=0)
    return float(np.diff(s.roots())) 

def fit_policy_function(self,PF):
        '''
        Fits the policy functions PF using the points xgrid using UnivariateSpline
        '''
        xgrid,S = self.xgrid,self.S
        
        Vf,cf,nf,xprimef = {},{},{},{}
        for s in range(S):
            PFvec = np.vstack(map(lambda x:PF(x,s),xgrid))
            Vf[s] = UnivariateSpline(xgrid,PFvec[:,0],s=0)
            cf[s] = UnivariateSpline(xgrid,PFvec[:,1],s=0,k=1)
            nf[s] = UnivariateSpline(xgrid,PFvec[:,2],s=0,k=1)
            for sprime in range(S):
                xprimef[s,sprime] = UnivariateSpline(xgrid,PFvec[:,3+sprime],s=0,k=1)
        
        return Vf,[cf,nf,xprimef] 

def process_trajectory(self, dxdt=True):
        """
        See https://docs.scipy.org/doc/scipy-0.15.1/reference/generated/scipy.interpolate.UnivariateSpline.html
        """
        self.trajectory_interp = []
        # Make sure s = 0.5 is ok.
        self.t_limits[:] = (np.min(self.input_history[:, 0]),
                            np.max(self.input_history[:, 0]))
        for i_dim in range(3):
            self.trajectory_interp.append(
                interpolate.UnivariateSpline(self.input_history[:, 0],
                                             self.input_history[:, i_dim + 1],
                                             k=1,
                                             s=0.,
                                             ext='raise'))
        if dxdt:
            self.trajectory_vel_interp = []
            for i_dim in range(3):
                self.trajectory_vel_interp.append(
                    self.trajectory_interp[i_dim].derivative()) 

def _residual(param, radial, profile, previous):
    """ `scipy.optimize.leastsq` residuals function.

        Evaluate the difference between a radial-scaled intensity profile
        and its adjacent "previous" angular slice.

    """

    radial_scaling, amplitude = param[0], param[1]

    newradial = radial * radial_scaling
    spline_prof = UnivariateSpline(newradial, profile, s=0, ext=3)
    newprof = spline_prof(radial) * amplitude

    # residual cf adjacent slice profile
    return newprof - previous 

def roc():
    '''
    Plots an ROC curve illustrating the performance of the classifier used in
    grain detection.
    '''
    n_classes = 2
    clf = get_model()
    (train_data, train_targets, test_data, test_targets) = Helper.unserialize("../Datasets/new_data_glcm_d1_a4_75_25.data")
    y_score = clf.decision_function(test_data)

    fpr, tpr, thresholds = metrics.roc_curve(test_targets, y_score)

    xnew = np.linspace(fpr.min(),fpr.max(),300)
    spl = UnivariateSpline(fpr,tpr)

    plt.figure()
    plt.plot(fpr, tpr, label='Exact ROC curve')
    plt.plot(xnew, spl(xnew), label='Smoothed ROC curve', color='red', linestyle='--')
    plt.xlim([0.0, 1.0])
    plt.ylim([0.0, 1.05])
    plt.xlabel('1 - Specificity (False Positive Rate)')
    plt.ylabel('Sensitivity (True Positive Rate)')
    plt.title('Receiver Operating Characteristic curve')
    plt.legend(loc="lower right")
    plt.show() 

def fit_inv_pdf(ladder_energies):
    """ Fit an interpolant to the inverse pdf of ladder energies
    Nasty hack to allow drawing energies from arbitrary distributions of ladder_energies

    Args:
      ladder_energies: array, output of ladder_numerical_err_hist

    Returns:
      interp_inv_pdf: inverse pdf interpolant
    """
    hist, bin_edges = np.histogram(ladder_energies, bins='auto')
    bin_mdpts = (np.diff(bin_edges) / 2) + bin_edges[:-1]
    # interpolation backward using slope bin_mdpts[1] - bin_mdpts[0]
    zero_interp_mdpt = -2 * bin_mdpts[0] + bin_mdpts[1]
    pdf = np.cumsum(hist) / np.sum(hist)
    # we add zero so that the interpolation is defined everywhere on [0, 1]
    pdf = np.concatenate([[0], pdf])
    bin_mdpts = np.concatenate([[zero_interp_mdpt], bin_mdpts])
    return UnivariateSpline(pdf, bin_mdpts, bbox=[0, 1], k=1) 

def on_iteration_end(self, iteration, status=None):
        y = status['unfolded']
        x = np.arange(len(y), dtype=float)
        if self.groups is None:
            spline = UnivariateSpline(x, y, k=self.degree, s=self.smooth)
            fitted_unfolded = spline(x)
        else:
            # Check that a group is specified for each cause bin
            if len(self.groups) != len(y):
                err_msg = ('Invalid groups array. There should be an entry '
                           'for each cause bin. However, got len(groups)={} '
                           'while there are {} cause bins.'.format(len(self.groups), len(y)))
                raise ValueError(err_msg)
            fitted_unfolded = np.empty(len(y))
            group_ids = np.unique(self.groups)
            for group in group_ids:
                group_mask = self.groups == group
                x_group = x[group_mask]
                y_group = y[group_mask]
                spline_group = UnivariateSpline(x_group, y_group, k=self.degree, s=self.smooth)
                fitted_unfolded_group = spline_group(x_group)
                fitted_unfolded[group_mask] = fitted_unfolded_group

        status['unfolded'] = fitted_unfolded 

def _integrate(self):
        """Walks through material list and computes the depth along the
        position (path). Computes the spline for the position-depth relation
        and determines the maximum depth for the material selection.

        Method does not need to be called by the user, instead the class
        calls it when necessary.
        """

        from scipy.interpolate import UnivariateSpline
        self.density_depth = None
        self.knots = [0.]
        self.X_int = [0.]

        for start, end, density, _ in self.mat_list:
            self.knots.append(end)
            self.X_int.append(density * (end - start) + self.X_int[-1])

        self._s_X2h = UnivariateSpline(self.X_int, self.knots, k=1, s=0.)
        self._s_h2X = UnivariateSpline(self.knots, self.X_int, k=1, s=0.)
        self._max_X = self.X_int[-1] 

def test_Zukauskas_tube_row_correction_refit():
    from scipy.interpolate import UnivariateSpline
    from ht.conv_tube_bank import Zukauskas_Czs_low_Re_staggered, Zukauskas_Czs_high_Re_staggered, Zukauskas_Czs_inline 
    # Commands used to obtain the fitted data:
    Zukauskas_Cz_Zs = [0.968219, 1.01968, 1.04164, 1.04441, 1.07539, 1.09332, 1.13914, 1.16636, 1.23636, 1.2394, 1.24505, 1.3125, 1.33358, 1.38554, 1.43141, 1.48282, 1.4876, 1.55352, 1.58004, 1.60466, 1.65726, 1.67493, 1.70188, 1.79682, 1.91823, 1.99323, 1.99665, 2.04002, 2.16306, 2.18556, 2.19045, 2.30691, 2.3086, 2.36006, 2.45272, 2.45413, 2.57543, 2.59826, 2.72341, 2.7451, 2.8896, 2.91482, 2.98759, 3.1572, 3.23203, 3.25334, 3.3511, 3.42295, 3.4499, 3.52072, 3.6168, 3.83565, 3.9076, 3.9826, 4.02939, 4.17411, 4.20042, 4.44242, 4.48937, 4.61023, 4.82811, 4.95071, 5.07038, 5.28825, 5.31232, 5.3621, 5.50606, 5.53014, 5.60405, 5.74801, 5.74807, 5.82181, 5.99012, 5.99017, 6.13636, 6.23207, 6.23212, 6.37826, 6.44983, 6.44988, 6.62015, 6.69183, 6.69188, 6.95807, 6.95812, 6.98312, 7.1767, 7.20001, 7.41772, 7.41848, 7.65967, 7.87743, 7.90156, 7.95003, 7.97416, 7.97476, 8.21606, 8.2166, 8.45795, 8.60365, 8.67571, 8.79712, 8.91809, 8.96597, 9.18368, 9.20824, 9.42551, 9.45013, 9.66741, 9.69197, 10.0786, 10.3208, 10.5623, 10.5626, 10.7803, 10.9737, 10.9978, 11.2398, 11.2399, 11.4574, 11.4575, 11.6993, 11.7478, 11.9653, 11.9896, 12.2072, 12.2315, 12.4491, 12.691, 12.7152, 12.9812, 13.2231, 13.2715, 13.465, 13.7068, 13.9246, 13.9487, 14.1905, 14.4324, 14.6743, 14.9161, 14.9887, 15.2305, 15.4724, 15.7142, 15.787, 15.811, 15.8835, 16.0046, 16.0287, 16.2465, 16.3673, 16.4883, 16.5124, 16.706, 16.7301, 16.9477, 16.9479, 17.1897, 17.2138, 17.4315, 17.6734, 17.9152, 17.9636, 18.2054, 18.2055, 18.4473, 18.6891, 18.9068, 18.931, 18.9793, 19.2212, 19.4631, 19.5599, 19.7049, 19.9467, 19.9952]
    low_Re_staggered_Cz = [0.828685, 0.831068, 0.832085, 0.832213, 0.833647, 0.834478, 0.836599, 0.83786, 0.8411, 0.841241, 0.841503, 0.845561, 0.84683, 0.849956, 0.852715, 0.855808, 0.856096, 0.859148, 0.860376, 0.861516, 0.863952, 0.864828, 0.866165, 0.870874, 0.876897, 0.880617, 0.880787, 0.882293, 0.886566, 0.887348, 0.887517, 0.89214, 0.892207, 0.894249, 0.897396, 0.897444, 0.901563, 0.902338, 0.906589, 0.907258, 0.911719, 0.912497, 0.914744, 0.91998, 0.92229, 0.922729, 0.92474, 0.926218, 0.926772, 0.928561, 0.930987, 0.936514, 0.938332, 0.940226, 0.940947, 0.943179, 0.943584, 0.946941, 0.947769, 0.9499, 0.95374, 0.955902, 0.957529, 0.960492, 0.96082, 0.961497, 0.962826, 0.963048, 0.96373, 0.965208, 0.965208, 0.965965, 0.967759, 0.96776, 0.969318, 0.969757, 0.969758, 0.970428, 0.970757, 0.970757, 0.971538, 0.972422, 0.972422, 0.975703, 0.975704, 0.976012, 0.978249, 0.978139, 0.977115, 0.977111, 0.977585, 0.978013, 0.97806, 0.978155, 0.978202, 0.978204, 0.97819, 0.97819, 0.979578, 0.980416, 0.980411, 0.980405, 0.981521, 0.981333, 0.980478, 0.980382, 0.981379, 0.981492, 0.981479, 0.981478, 0.982147, 0.982566, 0.982553, 0.982553, 0.98254, 0.981406, 0.98171, 0.98476, 0.984762, 0.98475, 0.98475, 0.984736, 0.984733, 0.985732, 0.985843, 0.986842, 0.986953, 0.985817, 0.986825, 0.986926, 0.986911, 0.987834, 0.988018, 0.988008, 0.987994, 0.991353, 0.991148, 0.98909, 0.9902, 0.990187, 0.991297, 0.991293, 0.991279, 0.991266, 0.992054, 0.992292, 0.99237, 0.992366, 0.992359, 0.992358, 0.993068, 0.993463, 0.993456, 0.993454, 0.994443, 0.994566, 0.994553, 0.994553, 0.99454, 0.994539, 0.996774, 0.99676, 0.996746, 0.996744, 0.99673, 0.99673, 0.997466, 0.998201, 0.998863, 0.998936, 0.99902, 0.999439, 0.999857, 1.00002, 1.00002, 1, 1]
    high_Re_staggered_Cz = [0.617923, 0.630522, 0.635897, 0.636344, 0.64134, 0.644232, 0.651621, 0.654452, 0.661728, 0.662045, 0.662632, 0.669643, 0.671835, 0.683767, 0.694302, 0.704706, 0.705673, 0.719014, 0.721221, 0.72327, 0.727649, 0.729119, 0.73359, 0.749337, 0.759443, 0.770509, 0.771014, 0.777413, 0.785006, 0.786394, 0.786756, 0.795376, 0.795545, 0.800697, 0.809975, 0.810062, 0.817547, 0.818955, 0.829084, 0.830839, 0.842534, 0.843935, 0.847977, 0.857398, 0.861555, 0.862739, 0.866619, 0.869471, 0.870563, 0.873432, 0.877325, 0.886614, 0.889668, 0.89251, 0.894282, 0.899765, 0.900781, 0.910119, 0.911931, 0.916595, 0.921077, 0.925619, 0.930052, 0.932064, 0.932286, 0.933053, 0.935273, 0.935644, 0.937165, 0.940127, 0.940128, 0.941835, 0.945731, 0.945731, 0.947081, 0.947964, 0.947965, 0.949465, 0.950199, 0.9502, 0.952562, 0.953557, 0.953558, 0.958036, 0.958037, 0.958267, 0.960054, 0.96027, 0.961381, 0.961388, 0.963615, 0.964614, 0.964725, 0.965472, 0.965844, 0.965847, 0.966954, 0.966957, 0.968064, 0.96956, 0.970299, 0.970762, 0.971224, 0.971406, 0.972518, 0.972516, 0.972504, 0.972617, 0.973614, 0.97368, 0.974715, 0.975264, 0.975811, 0.975814, 0.978048, 0.980033, 0.980281, 0.982515, 0.982515, 0.982502, 0.982503, 0.983612, 0.98361, 0.983597, 0.983709, 0.984707, 0.984819, 0.985817, 0.985804, 0.985896, 0.986911, 0.986898, 0.98712, 0.988008, 0.987994, 0.988994, 0.989104, 0.98909, 0.9902, 0.990187, 0.991297, 0.991293, 0.991279, 0.991266, 0.991252, 0.995742, 0.994903, 0.992366, 0.99573, 0.995729, 0.995716, 0.99571, 0.995703, 0.995826, 0.996814, 0.996813, 0.996801, 0.996801, 0.996787, 0.996786, 0.996774, 0.99676, 0.997682, 0.997867, 0.997854, 0.997854, 0.99784, 0.997826, 0.997814, 0.997813, 0.99781, 0.99892, 0.998907, 0.998901, 0.998893, 0.998879, 0.998877]
    inline_Cz = [0.658582, 0.681965, 0.69194, 0.6932, 0.700314, 0.704433, 0.710773, 0.714541, 0.724228, 0.724649, 0.725518, 0.735881, 0.738799, 0.74599, 0.751285, 0.75722, 0.757717, 0.76457, 0.767327, 0.776314, 0.781783, 0.783619, 0.786421, 0.794105, 0.803931, 0.81, 0.810227, 0.813093, 0.821227, 0.822615, 0.822917, 0.830103, 0.830207, 0.833383, 0.839101, 0.839188, 0.847046, 0.848103, 0.853898, 0.854902, 0.862547, 0.863881, 0.868371, 0.875104, 0.878568, 0.879555, 0.884081, 0.886933, 0.888003, 0.890813, 0.894236, 0.902032, 0.904114, 0.906285, 0.907639, 0.910812, 0.911389, 0.916696, 0.917725, 0.92029, 0.924912, 0.927513, 0.930052, 0.934534, 0.934906, 0.935673, 0.937893, 0.938227, 0.939252, 0.941249, 0.94125, 0.942957, 0.946853, 0.946854, 0.948204, 0.949088, 0.949088, 0.950588, 0.951322, 0.951323, 0.953685, 0.954679, 0.95468, 0.959159, 0.95916, 0.959274, 0.960163, 0.96027, 0.961381, 0.961388, 0.963615, 0.96585, 0.966222, 0.966969, 0.966968, 0.966968, 0.966954, 0.966957, 0.968064, 0.96956, 0.970299, 0.970762, 0.971224, 0.971406, 0.972518, 0.972516, 0.972504, 0.972617, 0.973614, 0.973737, 0.975675, 0.976888, 0.978099, 0.9781, 0.979191, 0.98016, 0.980281, 0.982515, 0.982515, 0.982502, 0.982503, 0.983612, 0.98361, 0.983597, 0.983709, 0.984707, 0.984819, 0.985817, 0.985804, 0.985896, 0.986911, 0.986898, 0.98712, 0.988008, 0.987994, 0.988994, 0.989104, 0.98909, 0.9902, 0.990187, 0.991297, 0.991293, 0.991279, 0.991266, 0.991252, 0.995742, 0.994903, 0.992366, 0.99573, 0.995729, 0.995716, 0.99571, 0.995703, 0.995826, 0.996814, 0.996813, 0.996801, 0.996801, 0.996787, 0.996786, 0.996774, 0.99676, 0.997682, 0.997867, 0.997854, 0.997854, 0.99784, 0.997826, 0.997814, 0.997813, 0.99781, 0.99892, 0.998907, 0.998901, 0.998893, 0.998879, 0.998877]
    
    # hand tuned smoothing
    Zukauskas_Cz_low_Re_staggered_obj = UnivariateSpline(Zukauskas_Cz_Zs, low_Re_staggered_Cz, s=0.0001) 
    Zukauskas_Cz_high_Re_staggered_obj = UnivariateSpline(Zukauskas_Cz_Zs, high_Re_staggered_Cz, s=0.0005)
    Zukauskas_Cz_inline_obj = UnivariateSpline(Zukauskas_Cz_Zs, inline_Cz, s=0.0005)
    
    Zukauskas_Czs_inline2 = np.round(Zukauskas_Cz_inline_obj(range(1, 20)), 4).tolist()
    assert_allclose(Zukauskas_Czs_inline, Zukauskas_Czs_inline2)

    Zukauskas_Czs_low_Re_staggered2 = np.round(Zukauskas_Cz_low_Re_staggered_obj(range(1, 20)), 4).tolist()
    assert_allclose(Zukauskas_Czs_low_Re_staggered, Zukauskas_Czs_low_Re_staggered2)

    Zukauskas_Czs_high_Re_staggered2 = np.round(Zukauskas_Cz_high_Re_staggered_obj(range(1, 20)), 4).tolist()
    assert_allclose(Zukauskas_Czs_high_Re_staggered, Zukauskas_Czs_high_Re_staggered2) 

def test_baffle_correction_Bell_fit():
    from ht.conv_tube_bank import Bell_baffle_configuration_tck
    # 125 us to create.
    spl = splrep(Bell_baffle_configuration_Fcs, Bell_baffle_configuration_Jcs, s=8e-5)
    [assert_allclose(i, j) for (i, j) in zip(spl, Bell_baffle_configuration_tck)]
    
    Bell_baffle_configuration_obj = UnivariateSpline(Bell_baffle_configuration_Fcs, 
                                                     Bell_baffle_configuration_Jcs, 
                                                     s=8e-5)
#    import matplotlib.pyplot as plt
#    plt.plot(Bell_baffle_configuration_Fcs, Bell_baffle_configuration_Jcs)
#    pts = np.linspace(0, 1, 5000)
#    plt.plot(pts, [Bell_baffle_configuration_obj(i) for i in pts])
#    plt.plot(pts, [0.55 + 0.72*i for i in pts]) # Serth and HEDH 3.3.6g misses the tip
#    plt.show()
# 

def smooth_elevation(df, smooth=4):
    if not present(df, Names.ELEVATION):
        log.debug(f'Smoothing {Names.RAW_ELEVATION} to get {Names.ELEVATION}')
        unique = df.loc[~df[Names.DISTANCE].isna() & ~df[Names.RAW_ELEVATION].isna(),
                        [Names.DISTANCE, Names.RAW_ELEVATION]].drop_duplicates(Names.DISTANCE)
        # the smoothing factor is from eyeballing results only.  maybe it should be a parameter.
        # it seems better to smooth along the route rather that smooth the terrain model since
        # 1 - we expect the route to be smoother than the terrain in general (roads / tracks)
        # 2 - smoothing the 2d terrain is difficult to control and can give spikes
        # 3 - we better handle errors from mismatches between terrain model and position
        #     (think hairpin bends going up a mountainside)
        # the main drawbacks are
        # 1 - speed on loading
        # 2 - no guarantee of consistency between routes (or even on the same routine retracing a path)
        spline = UnivariateSpline(unique[Names.DISTANCE], unique[Names.RAW_ELEVATION], s=len(unique) * smooth)
        df[Names.ELEVATION] = spline(df[Names.DISTANCE])
        df[Names.GRADE] = (spline.derivative()(df[Names.DISTANCE]) / 10)  # distance in km
        df[Names.GRADE] = df[Names.GRADE].rolling(5, center=True).median().ffill().bfill()
        # avoid extrapolation / interpolation
        df.loc[df[Names.RAW_ELEVATION].isna(), [Names.ELEVATION]] = None
    return df 

def _curvature_spline(skeletons, points_window=None, length=None):
    '''
    Calculate the curvature using univariate splines. This method is slower and can fail
    badly if the fit does not work, so I am only using it as testing
    '''

    def _spline_curvature(skel):
        if np.any(np.isnan(skel)):
            return np.full(skel.shape[0], np.nan)
        
        x = skel[:, 0]
        y = skel[:, 1]
        n = np.arange(x.size)

        fx = UnivariateSpline(n, x, k=5)
        fy = UnivariateSpline(n, y, k=5)

        x_d = fx.derivative(1)(n)
        x_dd = fx.derivative(2)(n)
        y_d = fy.derivative(1)(n)
        y_dd = fy.derivative(2)(n)

        curvature = _curvature_fun(x_d, y_d, x_dd, y_dd)
        return  curvature

    
    curvatures_fit = np.array([_spline_curvature(skel) for skel in skeletons])
    return curvatures_fit

#%% 

def moment_of_inertia(self):
        """
        Returns the moment of inertia of the layer [kg m^2]
        """
        moment = 0.0
        radii = self.radii
        density = self.evaluate(['density'])
        rhofunc = UnivariateSpline(radii, density)
        moment = np.abs(quad(lambda r: 8.0 / 3.0 * np.pi * rhofunc(r)
                             * r * r * r * r, radii[0], radii[-1])[0])
        return moment 

def mass(self):
        """
        Calculates the mass of the layer [kg]
        """
        mass = 0.0
        radii = self.radii
        density = self.evaluate(['density'])
        rhofunc = UnivariateSpline(radii, density)
        mass = np.abs(quad(lambda r: 4 * np.pi * rhofunc(r) *
                           r * r, radii[0], radii[-1])[0])
        return mass 

def test_dP_Kern_data():
    from ht.conv_tube_bank import Kern_f_Re_tck

    _Kern_dP_Res = np.array([9.9524, 11.0349, 12.0786, 13.0504, 14.0121, 15.0431, 16.1511, 17.1176, 17.9105, 18.9822,
        19.9879, 21.0484, 22.0217, 23.1893, 24.8973, 26.0495, 27.7862, 29.835, 31.8252, 33.9506, 35.9822, 38.3852,
        41.481, 43.9664, 47.2083, 50.6891, 54.0782, 58.0635, 63.5667, 68.2537, 74.247, 78.6957, 83.9573, 90.1511,
        95.5596, 102.613, 110.191, 116.806, 128.724, 137.345, 150.384, 161.484, 171.185, 185.031, 196.139, 210.639,
        230.653, 250.933, 281.996, 300.884, 329.472, 353.842, 384.968, 408.108, 444.008, 505.513, 560.821, 638.506,
        690.227, 741.254, 827.682, 918.205, 1018.63, 1122.76, 1213.62, 1320.38, 1417.94, 1522.93, 1667.69, 1838.11,
        2012.76, 2247.44, 2592.21, 2932.18, 3381.87, 3875.42, 4440.83, 5056.16, 5608.95, 6344.58, 7038.48, 8224.34,
        9123.83, 10121.7, 11598, 12701.4, 14090, 15938.5, 17452.9, 19112.6, 20929.3, 24614, 29324.6, 34044.8,
        37282.2, 42999.9, 50570.2, 55737.9, 59860.6, 65553, 70399.2, 78101.5, 84965.7, 96735.3, 110139, 122977,
        136431, 152339, 165740, 180319, 194904, 207981, 223357, 241440, 257621, 283946, 317042, 353996, 408315,
        452956, 519041, 590939, 668466, 751216, 827981, 894985, 1012440
    ])
    _Kern_dP_fs = 144.0 * np.array([0.0429177, 0.0382731, 0.0347901, 0.0316208, 0.0298653, 0.0276702, 0.0259671, 0.024523,
        0.0237582, 0.0224369, 0.0211881, 0.0202668, 0.0193847, 0.0184234, 0.0172894, 0.0166432, 0.0155182,
        0.0147509, 0.0138423, 0.0131572, 0.0124255, 0.0118105, 0.0110842, 0.0106028, 0.0100785, 0.00958019,
        0.0092235, 0.00871144, 0.00817649, 0.0077722, 0.00743616, 0.0071132, 0.00684836, 0.00655159, 0.00634789,
        0.00611185, 0.00592242, 0.00577517, 0.00552603, 0.00542355, 0.00522267, 0.00502847, 0.00493497, 0.00481301,
        0.00469334, 0.00460654, 0.00449314, 0.00438231, 0.00424799, 0.00416922, 0.00406658, 0.00401703, 0.00394314,
        0.0038947, 0.00382305, 0.00373007, 0.00368555, 0.00359592, 0.00357512, 0.003509, 0.00344515, 0.00338229,
        0.00332057, 0.00328077, 0.00322026, 0.00316102, 0.00308274, 0.00308446, 0.00302787, 0.00297247, 0.0028993,
        0.00284654, 0.00277759, 0.0027099, 0.00262738, 0.00256361, 0.00248541, 0.00244055, 0.00238072, 0.0023227,
        0.00228032, 0.00222531, 0.00218471, 0.00214484, 0.00206613, 0.00205439, 0.00200402, 0.00196775, 0.00191932,
        0.00189622, 0.00186143, 0.00180501, 0.0017393, 0.00170817, 0.00168761, 0.00163622, 0.00158663, 0.0015576,
        0.00153862, 0.0015201, 0.00149199, 0.00147418, 0.00142864, 0.00139389, 0.00136874, 0.00133524, 0.00131931,
        0.0012953, 0.00127147, 0.00124808, 0.00121724, 0.00121785, 0.00119533, 0.00118082, 0.00116638, 0.00114504,
        0.00111702, 0.00108969, 0.00107013, 0.00104389, 0.00101205, 0.000987437, 0.000969567, 0.000939849,
        0.000922653, 0.000905634, 0.000894962
    ])
#    # Used in preference over interp1d as saves 30% of execution time, and
#    # performs some marginally small amount of smoothing
#    # s=0.1 is chosen to have 9 knots, a reasonable amount.
#    Kern_f_Re = UnivariateSpline(_Kern_dP_Res, _Kern_dP_fs, s=0.1)
    tck = splrep(_Kern_dP_Res, _Kern_dP_fs, s=0.1)
    [assert_allclose(i, j) for i, j in zip(Kern_f_Re_tck, tck)] 

def test_Rac_Nusselt_Rayleigh_disk_fits():
    from fluids.optional import pychebfun
    from ht.conv_free_enclosed import insulated_disk_coeffs, uninsulated_disk_coeffs
    ratios = [0.4, 0.5, 0.7, 1.0, 1.4, 2.0, 3.0, 4.0, 6]
    Ras_uninsulated = [151200, 66600, 21300, 8010, 4350, 2540, 2010, 1880, 1708]
    Ras_insulated = [51800, 23800, 8420, 3770, 2650, 2260, 1900, 1830, 1708]
    
    uninsulated = UnivariateSpline(ratios, 1/np.log(Ras_uninsulated), k=1, s=0)
    insulated = UnivariateSpline(ratios, 1/np.log(Ras_insulated), k=1, s=0)
    
    N = 8
    insulated_fun = pychebfun.chebfun(insulated, domain=[ratios[0], ratios[-1]], N=N)
    uninsulated_fun = pychebfun.chebfun(uninsulated, domain=[ratios[0], ratios[-1]], N=N)
    

    insulated_coeffs = pychebfun.chebfun_to_poly(insulated_fun)
    uninsulated_coeffs = pychebfun.chebfun_to_poly(uninsulated_fun)
    
    assert_allclose(insulated_coeffs, insulated_disk_coeffs)
    assert_allclose(uninsulated_coeffs, uninsulated_disk_coeffs)
    
#    more_ratios = np.logspace(np.log10(ratios[0]), np.log10(ratios[-1]), 1000)
#    plt.semilogy(ratios, Ras_insulated)
#    plt.semilogy(ratios, Ras_uninsulated)
#    
#    plt.semilogy(more_ratios, np.exp(1/insulated_fun(np.array(more_ratios))), 'x')
#    plt.semilogy(more_ratios, np.exp(1/uninsulated_fun(np.array(more_ratios))), 'o')
#    plt.show() 

def test_SplineRegularizer_groups(example_dataset):
    degree = 3
    smooth = 20
    groups = np.empty_like(example_dataset.data)
    groups[:len(groups) // 2] = 0
    groups[len(groups) // 2:] = 1
    spline_reg = SplineRegularizer(degree=degree, smooth=smooth, groups=groups)
    unfolded_with_reg = iterative_unfold(data=example_dataset.data,
                                         data_err=example_dataset.data_err,
                                         response=example_dataset.response,
                                         response_err=example_dataset.response_err,
                                         efficiencies=example_dataset.efficiencies,
                                         efficiencies_err=example_dataset.efficiencies_err,
                                         return_iterations=True,
                                         callbacks=[spline_reg])

    unfolded_no_reg = iterative_unfold(data=example_dataset.data,
                                       data_err=example_dataset.data_err,
                                       response=example_dataset.response,
                                       response_err=example_dataset.response_err,
                                       efficiencies=example_dataset.efficiencies,
                                       efficiencies_err=example_dataset.efficiencies_err,
                                       return_iterations=True)
    # Manually regularize each group independently
    y_no_reg = unfolded_no_reg.iloc[0]['unfolded']
    x = np.arange(len(y_no_reg), dtype=float)
    fitted_unfolded_no_reg = np.empty(len(y_no_reg))
    group_ids = np.unique(groups)
    for group in group_ids:
        group_mask = groups == group
        x_group = x[group_mask]
        y_group = y_no_reg[group_mask]
        spline_group = UnivariateSpline(x_group, y_group, k=degree, s=smooth)
        fitted_unfolded_group = spline_group(x_group)
        fitted_unfolded_no_reg[group_mask] = fitted_unfolded_group

    np.testing.assert_allclose(unfolded_with_reg.iloc[0]['unfolded'],
                               fitted_unfolded_no_reg) 

def inv_transform_spline(self):
                '''Only for "nonstandard". Method to create inverse spline to discrete cumulative distribution function
                to allow drawing random variables.
                Warning: user should check that spline faithfully matches actual cdf.
                '''
                srt_param=np.sort(self.param)
                cdf = np.array(range(len(self.param)))/float(len(self.param))
                #create a spline
                self.spline2_cdf = UnivariateSpline(cdf,srt_param,k=5) 

def _solve(P_profile, T_profile, ref_pressure, mu_profile,
           planet_mass, planet_radius, star_radius,
           above_cloud_cond, T_star=None):
    assert(len(P_profile) == len(T_profile))

    # Ensure that the atmosphere is bound by making rough estimates of the
    # Hill radius and atmospheric height
    if T_star is None:
        T_star = Teff_sun

    R_hill = star_radius * \
        (T_star / T_profile[0])**2 * (planet_mass / (3 * M_sun))**(1.0 / 3)
    max_r_estimate = 1.0 / (1 / planet_radius + k_B * np.mean(T_profile) * np.log(
        P_profile[0] / P_profile[-1]) / (G * planet_mass * np.mean(mu_profile) * AMU))
    # The above equation is negative if the real answer is infinity

    if max_r_estimate < 0 or max_r_estimate > R_hill:
        raise AtmosphereError("Atmosphere unbound: height > hill radius")

    mu_interpolator = UnivariateSpline(np.log(P_profile), mu_profile, s=0)
    T_interpolator = UnivariateSpline(np.log(P_profile), T_profile, s=0)

    P_below = np.append(ref_pressure, P_profile[P_profile > ref_pressure])
    P_above = np.append(ref_pressure,
                        P_profile[P_profile <= ref_pressure][::-1])
    radii_below = _get_radii(np.log(P_below), planet_mass, planet_radius,
                             T_interpolator, mu_interpolator)
    radii_above = _get_radii(np.log(P_above), planet_mass, planet_radius,
                             T_interpolator, mu_interpolator)
    radii = np.append(radii_above.flatten()[1:][::-1],
                      radii_below.flatten()[1:])
    radii = radii[above_cloud_cond]
    dr = -np.diff(radii)

    return radii, dr 

def _interpolated_template(self, templateid):
        """Return an interpolator for the given template"""
        phase, y = self._get_template_by_id(templateid)

        # double-check that phase ranges from 0 to 1
        assert phase.min() >= 0
        assert phase.max() <= 1

        # at the start and end points, we need to add ~5 points to make sure
        # the spline & derivatives wrap appropriately
        phase = np.concatenate([phase[-5:] - 1, phase, phase[:5] + 1])
        y = np.concatenate([y[-5:], y, y[:5]])

        # Univariate spline allows for derivatives; use this!
        return UnivariateSpline(phase, y, s=0, k=5) 

def test_basic_template_model():
    template_id = 25

    try:
        templates = fetch_rrlyrae_templates()
    except(URLError, ConnectionError):
        raise SkipTest("No internet connection: "
                       "data download test skipped")

    phase, y = templates.get_template(templates.ids[template_id])
    model = UnivariateSpline(phase, y, s=0, k=5)

    theta = [17, 0.5, 0.3]
    period = 0.63
    rng = np.random.RandomState(0)
    t = rng.rand(20)
    mag = theta[0] + theta[1] * model((t / period - theta[2]) % 1)

    model = RRLyraeTemplateModeler('ugriz')
    model.fit(t, mag, 1)

    # check that the model matches what we expect
    assert_allclose(model._model(t, theta, period, template_id), mag)

    # check that the optimized model matches the input
    for use_gradient in [True, False]:
        theta_fit = model._optimize(period, template_id, use_gradient)
        assert_allclose(theta, theta_fit, rtol=1E-4)

    # check that the chi2 is near zero
    assert_allclose(model._chi2(theta_fit, period, template_id), 0,
                    atol=1E-8) 

def __init__(self, coefficient, frequency=None):
        if isinstance(coefficient, (int, float)):
            if frequency is not None and type(frequency) != float:
                coefficient = [coefficient for _ in range(len(frequency))]
            else:
                coefficient = [coefficient]

        self.coefficient = coefficient
        self.frequency = frequency

        if len(self.coefficient) > 1:
            try:
                with warnings.catch_warnings():
                    warnings.simplefilter("ignore")
                    self.interpolated = interpolate.UnivariateSpline(
                        self.frequency, self.coefficient
                    )
            #  dfitpack.error is not exposed by scipy
            #  so a bare except is used
            except:
                try:
                    if len(self.frequency) in (2, 3):
                        self.interpolated = interpolate.interp1d(
                            self.frequency,
                            self.coefficient,
                            kind=len(self.frequency) - 1,
                            fill_value="extrapolate",
                        )
                except:
                    raise ValueError(
                        "Arguments (coefficients and frequency)"
                        " must have the same dimension"
                    )
        else:
            self.interpolated = lambda x: np.array(self.coefficient[0]) 

def set_virtual(self, motors_maximum_powers):
        from scipy.interpolate import UnivariateSpline as Spline
        p, pb = self.hev_power_model.delta_ice_power(motors_maximum_powers)[:-1]
        pb, i = np.unique(-pb, return_index=True)
        self._battery_power = Spline(pb, np.abs(p.ravel()[i]), k=1, s=0)

    # noinspection PyUnresolvedReferences 

def bellman_operator(pgrid, c, f0, f1, L0, L1, J):
    """
    Evaluates the value function for a given continuation value
    function; that is, evaluates

        J(p) = min((1 - p) L0, p L1, c + E J(p'))

    Uses linear interpolation between points.
    """
    m = np.size(pgrid)
    assert m == np.size(J)
    
    J_out = np.zeros(m)
    J_interp = interp.UnivariateSpline(pgrid, J, k=1, ext=0)

    for (p_ind, p) in enumerate(pgrid):
        # Payoff of choosing model 0
        p_c_0 = expect_loss_choose_0(p, L0)
        p_c_1 = expect_loss_choose_1(p, L1)
        p_con = expect_loss_cont(p, c, f0, f1, J_interp)
        
        J_out[p_ind] = min(p_c_0, p_c_1, p_con)

    return J_out


#  == Now run at given parameters == #

#  First set up distributions