def METHOD_NAME(records, time_range, seconds_range, config, to_pe,
                   max_samples=DEFAULT_MAX_SAMPLES,
                   ignore_max_sample_warning=False):
    """Return (wv_matrix, times, pms)
      - wv_matrix: (n_samples, n_pmt) array with per-PMT waveform intensity in PE/ns
      - times: time labels in seconds (corr. to rows)
      - pmts: PMT numbers (corr. to columns)
    Both times and pmts have one extra element.
    :param max_samples: Maximum number of time samples. If window and dt
    conspire to exceed this, waveforms will be downsampled.
    :param ignore_max_sample_warning: If True, suppress warning when this happens.
    Example:
        wvm, ts, ys = st.records_matrix(run_id, seconds_range=(1., 1.00001))
        plt.pcolormesh(ts, ys, wvm.T,
                       norm=matplotlib.colors.LogNorm())
        plt.colorbar(label='Intensity [PE / ns]')
    """
    if len(records):
        dt = records[0]['dt']
        samples_per_record = len(records[0]['data'])
    else:
        # Defaults here do not matter, nothing will be plotted anyway
        dt = 10, 110
    record_duration = samples_per_record * dt
    window = time_range[1] - time_range[0]
    if window / dt > max_samples:
        with np.errstate(divide='ignore', invalid='ignore'):
            # Downsample. New dt must be
            #  a) multiple of old dt
            dts = np.arange(0, record_duration + dt, dt).astype(np.int64)
            #  b) divisor of record duration
            dts = dts[record_duration / dts % 1 == 0]
            #  c) total samples < max_samples
            dts = dts[window / dts < max_samples]
            if len(dts):
                # Pick lowest dt that satisfies criteria
                dt = dts.min()
            else:
                # Records will be downsampled to single points
                dt = max(record_duration, window // max_samples)
        if not ignore_max_sample_warning:
            warnings.warn(f"Matrix would exceed max_samples {max_samples}, "
                          f"downsampling to dt = {dt} ns.")
    wvm = _records_to_matrix(
        records,
        t0=time_range[0],
        n_channels=config['n_tpc_pmts'],
        dt=dt,
        window=window)
    wvm = wvm.astype(np.float32) * to_pe.reshape(1, -1) / dt
    # Note + 1, so data for sample 0 will range from 0-1 in plot
    ts = (np.arange(wvm.shape[0] + 1) * dt / int(1e9) + seconds_range[0])
    ys = np.arange(wvm.shape[1] + 1)
    return wvm, ts, ys