<SYSTEM_TASK:>
Solve the following problem using Python, implementing the functions described below, one line at a time
<END_TASK>
<USER_TASK:>
Description:
 def bandstats(filenames=None, num_sample_points=3, temperature=None, degeneracy_tol=1e-4, parabolic=True):
 """Calculate the effective masses of the bands of a semiconductor. Args: filenames (:obj:`str` or :obj:`list`, optional):
 Path to vasprun.xml or vasprun.xml.gz file. If no filenames are provided, the code will search for vasprun.xml or vasprun.xml.gz files in folders named 'split-0*'. Failing that, the code will look for a vasprun in the current directory. If a :obj:`list` of vasprun files is provided, these will be combined into a single band structure. num_sample_points (:obj:`int`, optional):
 Number of k-points to sample when fitting the effective masses. temperature (:obj:`int`, optional):
 Find band edges within kB * T of the valence band maximum and conduction band minimum. Not currently implemented. degeneracy_tol (:obj:`float`, optional):
 Tolerance for determining the degeneracy of the valence band maximum and conduction band minimum. parabolic (:obj:`bool`, optional):
 Use a parabolic fit of the band edges. If ``False`` then nonparabolic fitting will be attempted. Defaults to ``True``. Returns: dict: The hole and electron effective masses. Formatted as a :obj:`dict` with keys: ``'hole_data'`` and ``'electron_data'``. The data is a :obj:`list` of :obj:`dict` with the keys: 'effective_mass' (:obj:`float`) The effective mass in units of electron rest mass, :math:`m_0`. 'energies' (:obj:`numpy.ndarray`) Band eigenvalues in eV. 'distances' (:obj:`numpy.ndarray`) Distances of the k-points in reciprocal space. 'band_id' (:obj:`int`) The index of the band, 'spin' (:obj:`~pymatgen.electronic_structure.core.Spin`) The spin channel 'start_kpoint' (:obj:`int`) The index of the k-point at which the band extrema occurs 'end_kpoint' (:obj:`int`) """

    if not filenames:
        filenames = find_vasprun_files()
    elif isinstance(filenames, str):
        filenames = [filenames]

    bandstructures = []
    for vr_file in filenames:
        vr = BSVasprun(vr_file, parse_projected_eigen=False)
        bs = vr.get_band_structure(line_mode=True)
        bandstructures.append(bs)
    bs = get_reconstructed_band_structure(bandstructures)

    if bs.is_metal():
        logging.error('ERROR: System is metallic!')
        sys.exit()

    _log_band_gap_information(bs)

    vbm_data = bs.get_vbm()
    cbm_data = bs.get_cbm()

    logging.info('\nValence band maximum:')
    _log_band_edge_information(bs, vbm_data)

    logging.info('\nConduction band minimum:')
    _log_band_edge_information(bs, cbm_data)

    if parabolic:
        logging.info('\nUsing parabolic fitting of the band edges')
    else:
        logging.info('\nUsing nonparabolic fitting of the band edges')

    if temperature:
        logging.error('ERROR: This feature is not yet supported!')

    else:
        # Work out where the hole and electron band edges are.
        # Fortunately, pymatgen does this for us. Points at which to calculate
        # the effective mass are identified as a tuple of:
        # (spin, band_index, kpoint_index)
        hole_extrema = []
        for spin, bands in vbm_data['band_index'].items():
            hole_extrema.extend([(spin, band, kpoint) for band in bands
                                 for kpoint in vbm_data['kpoint_index']])

        elec_extrema = []
        for spin, bands in cbm_data['band_index'].items():
            elec_extrema.extend([(spin, band, kpoint) for band in bands
                                 for kpoint in cbm_data['kpoint_index']])

        # extract the data we need for fitting from the band structure
        hole_data = []
        for extrema in hole_extrema:
            hole_data.extend(get_fitting_data(bs, *extrema,
                             num_sample_points=num_sample_points))

        elec_data = []
        for extrema in elec_extrema:
            elec_data.extend(get_fitting_data(bs, *extrema,
                             num_sample_points=num_sample_points))

    # calculate the effective masses and log the information
    logging.info('\nHole effective masses:')
    for data in hole_data:
        eff_mass = fit_effective_mass(data['distances'], data['energies'],
                                      parabolic=parabolic)
        data['effective_mass'] = eff_mass
        _log_effective_mass_data(data, bs.is_spin_polarized, mass_type='m_h')

    logging.info('\nElectron effective masses:')
    for data in elec_data:
        eff_mass = fit_effective_mass(data['distances'], data['energies'],
                                      parabolic=parabolic)
        data['effective_mass'] = eff_mass
        _log_effective_mass_data(data, bs.is_spin_polarized)

    return {'hole_data': hole_data, 'electron_data': elec_data}