def estimateDerivatives_fn(xS, startInds, numDtStepsForDeriv, pointWtsForThisVar, dt):
    """
    For a vector, calculate the target derivative (as in dx/dt = rise/run) that we hope to match
    when we regress on the library functionals. The target derivative is calculated using the
    variables' time-series. Also return the weights for each element of the target rise, for use
    in the regression.
    NOTE: We require values in xS before the first startInd and beyond the last start ind

    Parameters
    ----------
    xS : np.array of floats, (column vector same length as timepoints eg tTrain), the
    pre-processed time-series of one variable.
    startInds : np.array (vector of ints). Indices of timepoints.
    numDtStepsForDeriv : scalar int
    pointWtsForThisVar :  np.array (vector of floats, same length as timepoints eg tTrain).
    dt : scalar float. The length of the simulation timestep

    Returns
    -------
    derivEstimate : np.array (vector of floats). size = startInds.shape
    weights : np.array (vector of floats). size = startInds.shape

    """
    timeStep = dt*numDtStepsForDeriv  # float (ie a measure of time, not an index of timepoints)

    # 4th order approx to derivative:
    m2Inds = startInds - 2*numDtStepsForDeriv  # indices
    m1Inds = startInds - 1*numDtStepsForDeriv
    p2Inds = startInds + 2*numDtStepsForDeriv
    p1Inds = startInds + 1*numDtStepsForDeriv
    derivEstimates = \
        (xS[m2Inds] - 8*xS[m1Inds] + 8*xS[p1Inds] - xS[p2Inds]) / (12*timeStep)
    # The weight for each deriv estimate combines the timepoint values used:
    weights = \
        pointWtsForThisVar[m2Inds] + 8 * pointWtsForThisVar[m1Inds] + \
            8 * pointWtsForThisVar[p1Inds] + pointWtsForThisVar[p2Inds]
    weights = weights / sum(weights)

    return derivEstimates, weights