ts - An AstonSeries that's being fitted with peaks
    f - The functional form of the peaks (e.g. gaussian)
    rts - peak maxima to fit; each number corresponds to one peak

def guess_initc(ts, f, rts=[]):
    """
    ts - An AstonSeries that's being fitted with peaks
    f - The functional form of the peaks (e.g. gaussian)
    rts - peak maxima to fit; each number corresponds to one peak
    """
    def find_side(y, loc=None):
        if loc is None:
            loc = y.argmax()
        ddy = np.diff(np.diff(y))
        lft_loc, rgt_loc = loc - 2, loc + 1
        while rgt_loc >= 0 and rgt_loc < len(ddy):
            if ddy[rgt_loc] < ddy[rgt_loc - 1]:
                break
            rgt_loc += 1
        while lft_loc >= 0 and lft_loc < len(ddy):
            if ddy[lft_loc] < ddy[lft_loc + 1]:
                break
            lft_loc -= 1
        return lft_loc + 1, rgt_loc + 1

    # weight_mom = lambda m, a, w: \
    #   np.sum(w * (a - np.sum(w * a) / np.sum(w)) ** m) / np.sum(w)
    # sig = np.sqrt(weight_mom(2, ts.index, ts.values))  # sigma
    # peak_params['s'] = weight_mom(3, ts.index, ts.values) / sig ** 3
    # peak_params['e'] = weight_mom(4, ts.index, ts.values) / sig ** 4 - 3
    # TODO: better method of calculation of these?
    all_params = []
    for rt in rts:
        peak_params = {'x': rt}  # ts.index[ts.values.argmax()]
        top_idx = np.abs(ts.index - rt).argmin()
        side_idx = find_side(ts.values, top_idx)
        peak_params['h'] = ts.values[top_idx]
        # - min(ts.y[side_idx[0]], ts.y[side_idx[1]])
        peak_params['w'] = ts.index[side_idx[1]] - ts.index[side_idx[0]]
        peak_params['s'] = 1.1
        peak_params['e'] = 1.
        peak_params['a'] = 1.
        all_params.append(peak_params)
    return all_params