class BaseFeature(object):
    """
    A base class for all the features.

    Parameters
    ----------
    input_type : str, list of str, or callable, default='list'
        Specifies the format the input values will be (must be one of 'list',
        'filename', a list of strings, or a callable). If it is a list of
        strings, the strings tell the order of (and if they are included) the
        different molecule attributes (coords, elements, numbers,
        connections). If a callable is given, then it is assumed to return a
        LazyValues object.

    n_jobs : int, default=1
        Specifies the number of processes to create when generating the
        features. Positive numbers specify a specifc amount, and numbers less
        than 1 will use the number of cores the computer has.
    """
    def __init__(self, input_type='list', n_jobs=1):
        self.input_type = input_type
        self.n_jobs = n_jobs

    def _get_param_strings(self):
        argspec = inspect.getargspec(type(self).__init__)
        # Delete the only non-keyword argument
        args = [x for x in argspec.args if x != "self"]
        values = [getattr(self, x) for x in args]
        return ["%s=%r" % (x, y) for x, y in zip(args, values)]

    def __repr__(self):
        name = type(self).__name__
        params = self._get_param_strings()
        return "%s(%s)" % (name, ', '.join(params))

    def set_params(self, **kwargs):
        """
        Set the feature parameter values.

        Parameters
        ----------
            kwargs : kwargs
                Key value pairs to set for the feature parameters. Keys that
                are not valid parameters will be ignored.
        """
        for key, value in kwargs.items():
            try:
                getattr(self, key)
                setattr(self, key, value)
            except AttributeError:
                continue

    def get_params(self):
        """
        Get a dictonary of all the feature parameters.

        Returns
        -------
            params : dict
                A dictonary of all the feature parameters.
        """
        argspec = inspect.getargspec(type(self).__init__)
        # Delete the only non-keyword argument
        args = [x for x in argspec.args if x != "self"]
        values = [getattr(self, x) for x in args]
        return {key: value for key, value in zip(args, values)}

    def slugify(self):
        """
        Convert an instance to a simple string.

        Returns
        -------
        string : str
            The slug string

        """
        name = type(self).__name__
        # Skip the first two parameters
        params = self._get_param_strings()[2:]
        string = '__'.join([name] + params).replace("'", '')
        return string

    def convert_input(self, X):
        """
        Convert the input (as specified in self.input_type) to a usable form.

        Parameters
        ----------
        X : list or string (depends on the instance value of input_type)
            An object that stores the data for a single molecule. See the
            Notes for more details.

        Returns
        -------
        values : Object
            An object that allows the lazy evaluation of different properties

        Raises
        ------
        ValueError
            If the input_type given is not allowed.

        Notes
        -----
        If input_type is 'list', then it must be an iterable of (elements,
        coodinates pairs) for each molecule. Where the elements are an
        iterable of the form (ele1, ele2, ..., elen) and coordinates are an
        iterable of the form [(x1, y1, z1), (x2, y2, z2), ..., (xn, yn,
        zn)]. This allows allows for connections to be included. This is a
        dictionary where the keys are the indices of the atoms and the
        values are dictonaries with the key being another index and the
        value is the bond order (one of '1', 'Ar', '2', or '3').
        Example for methane::

            {
                0: {1: "1", 2: "1", 3: "1", 4: "1"},
                1: {0: "1"},
                2: {0: "1"},
                3: {0: "1"},
                4: {0: "1"},
            }

        If input_type is 'filename', then it must be an iterable of
        paths/filenames for each molecule. Currently, the supported formats
        are: xyz, mol2, and a simple xyz format (.out).

        If input_type is a list, then they will be treated as labels to
        each of the arguments passed in via a tuple. For example,
        input_type="list" can be reproduced with ["elements", "coords"]
        or ["elements", "coords", "connections"].

        If input_type is a callable, then it is assumed that the callable
        returns a LazyValues object.
        """
        connections = None
        if self.input_type == "list":
            try:
                first, coords = X
            except ValueError:
                first, coords, connections = X
            if len(first) and isinstance(first[0], str):
                values = LazyValues(elements=first, coords=coords,
                                    connections=connections)
            else:
                values = LazyValues(numbers=first, coords=coords,
                                    connections=connections)
        elif self.input_type == "filename":
            values = read_file_data(X)
        elif type(self.input_type) in (list, tuple):
            d = {x: y for x, y in zip(self.input_type, X)}
            values = LazyValues(**d)
        elif callable(self.input_type):
            values = self.input_type(X)
        else:
            raise ValueError("The input_type '%s' is not allowed." %
                             self.input_type)
        return values

    def get_labels(self):
        """
        Get the labels for the features in the transformer

        Returns
        -------
        values : tuple
            All of the labels of the resulting features.
        """
        if self.LABELS is None:
            return tuple()

        values = []
        for x in self.LABELS:
            try:
                func_name, data_name = x
                func = getattr(self, func_name)
                if data_name is None:
                    temp = func()
                else:
                    data = getattr(self, data_name)
                    temp = func(data)
            except (TypeError, ValueError) as e:
                if len(x) == 2:
                    # Error in calling the function
                    raise e
                temp = getattr(self, x)
            values.append(tuple(temp))
        return sum(values, tuple())

    def check_fit(self):
        """
        Check if the transformer has been fit

        Raises
        ------
        ValueError
            The transformer has not been fit.
        """
        if self.ATTRIBUTES is None:
            return

        msg = "This %s instance is not fitted yet. Call 'fit' first."
        for key in self.ATTRIBUTES:
            if getattr(self, key) is None:
                raise ValueError(msg % type(self).__name__)

    @classmethod
    def get_citation(self):
        try:
            docs = self.__doc__
            idx = docs.index("References")
            values = [x.strip() for x in docs[idx:].split('\n')[2:]]
            new_values = [[]]
            for value in values:
                if "----" in value:
                    new_values.pop()
                    break
                elif not value:
                    new_values.append([])
                else:
                    new_values[-1].append(value)
            strings = [' '.join(x) for x in new_values if x]
            return '\n'.join(strings)
        except ValueError:
            return "MolML https://github.com/crcollins/molml"

    def map(self, f, seq):
        """
        Parallel implementation of map.

        Parameters
        ----------
        f : callable
            A function to map to all the values in 'seq'

        seq : iterable
            An iterable of values to process with 'f'

        Returns
        -------
        results : list, shape=[len(seq)]
            The evaluated values
        """
        if self.n_jobs < 1:
            n_jobs = multiprocessing.cpu_count()
        elif self.n_jobs == 1:
            return list(map(f, seq))
        else:
            n_jobs = self.n_jobs

        pool = Pool(n_jobs)
        results = list(pool.map(f, seq))
        # Closing/joining is not really allowed because pathos sees pools as
        # lasting for the duration of the program.
        return results

    def reduce(self, f, seq):
        """
        Parallel implementation of reduce.

        This changes the problem from being O(n) steps to O(lg n)

        Parameters
        ----------
        f : callable
            A function to use to reduce the values of 'seq'

        seq : iterable
            An iterable of values to process

        Returns
        -------
        results : object
            A single reduced object based on 'seq' and 'f'
        """
        if self.n_jobs == 1:
            return reduce(f, seq)

        while len(seq) > 1:
            pairs = [(f, x, y) for x, y in zip(seq[::2], seq[1::2])]
            temp_seq = self.map(_func_star, pairs)
            # If the sequence length is odd add the last element on
            # This is because it will not get included with the zip
            if len(seq) % 2:
                temp_seq.append(seq[-1])
            seq = temp_seq
        return seq[0]

    def fit(self, X, y=None):
        """
        Fit the model.

        Parameters
        ----------
        X : list, shape=(n_samples, )
            A list of objects to use to fit.

        Returns
        -------
        self : object
            Returns the instance itself.
        """
        raise NotImplementedError

    def _para_transform(self, X):
        """
        A single instance of the transform procedure.

        This is formulated in a way that the transformations can be done
        completely parallel with map.

        Parameters
        ----------
        X : object
            An object to use for the transform

        Returns
        -------
        value : array-like
            The features extracted from the molecule
        """
        raise NotImplementedError

    def transform(self, X, y=None):
        """
        Framework for a potentially parallel transform.

        Parameters
        ----------
        X : list, shape=(n_samples, )
            A list of objects to use to transform

        Returns
        -------
        array : array, shape=(n_samples, n_features)
            The transformed features
        """
        self.check_fit()
        results = self.map(self._para_transform, X)
        return numpy.array(results)

    def fit_transform(self, X, y=None):
        """
        A naive default implementation of fitting and transforming.

        Parameters
        ----------
        X : list, shape=(n_samples, )
            A list of objects to use to fit and then transform

        Returns
        -------
        array : array, shape=(n_samples, n_features)
            The transformed features
        """
        self.fit(X, y)
        return self.transform(X, y)

    def to_json(self):
        """
        Return model data as a json compatible dict

        This will recursively convert other transformer objects as well.

        Returns
        -------
        data : dict
            The json data
        """
        attributes = {}
        if self.ATTRIBUTES is not None:
            attributes = {key: getattr(self, key) for key in self.ATTRIBUTES}

        full_name = self.__module__ + '.' + self.__class__.__name__
        params = {}
        for key, value in self.get_params().items():
            try:
                params[key] = value.to_json()
            except AttributeError:
                params[key] = value

        data = {
                "transformer": full_name,
                "parameters": params,
                "attributes": attributes,
        }
        return data

    def save_json(self, f):
        """
        Save the model data in a json file

        Parameters
        ----------
        f : str or file descriptor
            The path to save the data or a file descriptor to save it to.
        """
        data = self.to_json()
        try:
            json.dump(data, f)
        except AttributeError:
            with open(f, 'w') as out_file:
                json.dump(data, out_file) class SetMergeMixin(object):
    """
    A simple mixin that will merge sets.

    This mixin replaces all the duplicate code that just merges sets when
    doing the parallel fits. For this to work, it requires that the subclasses
    define `ATTRIBUTES`.
    """
    def fit(self, X, y=None):
        """
        Fit the model.

        Parameters
        ----------
        X : list, shape=(n_samples, )
            A list of objects to use to fit.

        Returns
        -------
        self : object
            Returns the instance itself.
        """
        res = self.map(self._para_fit, X)
        if len(self.ATTRIBUTES) > 1:
            temp = self.reduce(lambda x, y: tuple(set(xx) | set(yy)
                                                  for xx, yy in zip(x, y)),
                               res)
            for attr, vals in zip(self.ATTRIBUTES, temp):
                setattr(self, attr, tuple(sorted(set(vals))))
        else:
            vals = self.reduce(lambda x, y: set(x) | set(y), res)
            setattr(self, self.ATTRIBUTES[0], tuple(sorted(set(vals))))
        return self class EncodedFeature(BaseFeature):
    """
    This is a generalized class to handle all kinds of encoding feature
    representations. These approaches seem to be a fairly general way of
    making lists of scalar values more effective to use in machine learning
    models. Essentially, it can be viewed as kernel smoothed histograms over
    the values of interest.

    Parameters
    ----------
    input_type : string, default='list'
        Specifies the format the input values will be (must be one of 'list'
        or 'filename').

    n_jobs : int, default=1
        Specifies the number of processes to create when generating the
        features. Positive numbers specify a specifc amount, and numbers less
        than 1 will use the number of cores the computer has.

    segments : int, default=100
        The number of bins/segments to use when generating the histogram.
        Empirically, it has been found that values beyond 50-100 have little
        benefit.

    smoothing : string or callable, default='norm'
        A string or callable to use to smooth the histogram values. If a
        callable is given, it must take just a single argument that is a float
        (or vector of floats). For a list of supported default functions look
        at SMOOTHING_FUNCTIONS.

    start : float, default=0.2
        The starting point for the histgram sampling in angstroms.

    end : float, default=6.0
        The ending point for the histogram sampling in angstroms.

    slope : float, default=20.
        A parameter to tune the smoothing values. This is applied as a
        multiplication before calling the smoothing function.

    spacing : string or callable, default='linear'
        The histogram interval spacing type. Must be one of ("linear",
        "inverse", or "log"). Linear spacing is normal spacing. Inverse takes
        and evaluates the distances as 1/r and the start and end points are
        1/x. For log spacing, the distances are evaluated as numpy.log(r)
        and the start and end points are numpy.log(x). If the value is
        callable, then it should take a float or vector of floats and return
        a similar mapping like the other methods.

    References
    ----------
    Collins, C.; Gordon, G.; von Lilienfeld, O. A.; Yaron, D. Constant Size
    Molecular Descriptors For Use With Machine Learning. arXiv:1701.06649
    """
    def __init__(self, input_type='list', n_jobs=1, segments=100,
                 smoothing='norm', slope=20., start=0.2, end=6.,
                 spacing='linear'):
        super(EncodedFeature, self).__init__(input_type=input_type,
                                             n_jobs=n_jobs)
        self.segments = segments
        self.smoothing = smoothing
        self.slope = slope
        self.start = start
        self.end = end
        self.spacing = spacing

    def _get_theta_info(self):
        theta_func = get_spacing_function(self.spacing)
        theta = numpy.linspace(theta_func(self.start), theta_func(self.end),
                               self.segments)
        return theta, theta_func

    def encode_values(self, iterator, lengths, saved_lengths=0):
        '''
        Encodes an iterable of values into a uniform length array. These
        values can then be indexed to allow binning them in different sections
        of the array. After the values are processed, the array can by
        flattened down to a desired number of axes.

        Parameters
        ----------
            iterator : iterable
                The collection of values to encode. Each item in the iterable
                must contain values for (idx, value, scaling). Where idx is a
                tuple of integer values indicating which encoding bucket the
                values go in, value is the value to encode, and scaling is a
                factor that gets multiplied by the final encoded subvector
                before getting added to the total (This is mostly used to mask
                values and scale their influence with distance. If idx is None,
                then the value will be skipped.

            lengths : tuple of ints
                The number of encoding axes to create. In terms of
                EncodedBonds, this would be the number of element pairs.

            saved_lengths : ints
                The number of axis components to retain. The order that they
                get saved is the same order that is given in lengths. For
                example, when doing atom encodings, this should be 1 to retain
                the atom axis.

        Returns
        -------
            vector : array
                The final concatenated vector of all the subvectors. This will
                have a shape of lengths[:saved_lengths]
                + product(lengths[saved_lengths:]) * segments).
        '''
        smoothing_func = get_smoothing_function(self.smoothing)
        vector = numpy.zeros(tuple(lengths) + (self.segments, ))
        theta, theta_func = self._get_theta_info()

        for idxs, value, scaling in iterator:
            if idxs is None:
                continue
            diff = theta - theta_func(value)
            value = smoothing_func(diff, self.slope)
            vector[tuple(idxs)] += value * scaling

        reshape = tuple(lengths)[:saved_lengths] + (-1, )
        return vector.reshape(*reshape)

    def get_group_order(self, groups):
        """
        Parameters
        ----------
        groups : list
            A list of all the groups.

        Returns
        -------
        value_order : list
            A list of all groups in order.
        """
        return groups

    def get_encoded_labels(self, groups):
        """
        Parameters
        ----------
        groups : list
            A list of all the groups.

        Returns
        -------
        labels : list
            A list of all the feature labels.
        """
        theta, theta_func = self._get_theta_info()
        labels = []
        for group in self.get_group_order(groups):
            name = '-'.join(group)
            for x in theta:
                labels.append('%s_%s' % (name, round(x, 5)))
        return labels class FormMixin(object):
    """
    A simple mixin for handling form transformations

    This mixin handles all how index mapping is done when going from higher
    dimensional attributes to lower dimensional ones. By default, this mixin
    uses the first value in ATTRIBUTES as the basis for the index mapping.

    use_comb_idxs : bool, default=False
        Whether or not to use all combinations of indices when doing the
        subselection. If this is false, a  middle out scheme will be used.
    """
    def __init__(self, form=1, add_unknown=False, use_comb_idxs=False,
                 *args, **kwargs):
        super(FormMixin, self).__init__(*args, **kwargs)
        self.form = form
        self.add_unknown = add_unknown
        self.use_comb_idxs = use_comb_idxs
        self._idx_map = None

    def get_idx_map(self):
        """
        Lazily load the idx_map.

        Returns
        -------
        idx_map : IndexMap
            The IndexMap object for this form and add_unknown setting.
        """
        values = getattr(self, self.ATTRIBUTES[0])
        if self._idx_map is None or not self._idx_map.is_valid(values):
            self._idx_map = IndexMap(values, self.form, self.add_unknown,
                                     self.use_comb_idxs)
        return self._idx_map

    def get_group_order(self, groups):
        """
        Parameters
        ----------
        groups : list
            A list of all the groups. This is ignored.

        Returns
        -------
        value_order : list
            A list of all groups in order.
        """
        return self.get_idx_map().get_value_order()

    def transform(self, X, y=None):
        """
        Framework for a potentially parallel transform.

        Parameters
        ----------
        X : list, shape=(n_samples, )
            A list of objects to use to transform

        Returns
        -------
        array : array, shape=(n_samples, n_features)
            The transformed features
        """
        self.check_fit()
        # This is to ensure that the idx_map exists before potentially
        # duplicating the process for parallelization.
        self.get_idx_map()
        results = self.map(self._para_transform, X)
        return numpy.array(results) def get_depth_threshold_mask_connections(connections, min_depth=0,
                                         max_depth=numpy.inf):
    """
    Get the depth threshold mask from connections.

    Parameters
    ----------
    connections : dict, index->list of indices
        A dictionary that contains lists of all connected atoms.


    min_depth : int, default=0
        The minimum depth to allow in the masking

    max_depth : int, default=numpy.inf
        The maximum depth to allow in the masking

    Returns
    -------
    mask : numpy.array, shape=(len(connections), len(connections))
        A mask of all the atoms that are less than or equal to `max_depth`
        away.
    """
    if max_depth < 1:
        max_depth = numpy.inf
    if not min_depth and max_depth is numpy.inf:
        V = len(connections)
        return numpy.ones((V, V)).astype(bool)
    dist = get_graph_distance(connections)
    return (min_depth <= dist) & (dist <= max_depth) def get_coulomb_matrix(numbers, coords, alpha=1, use_decay=False):
    r"""
    Return the coulomb matrix for the given coords and numbers.

    .. math::

        C_{ij} = \begin{cases}
            \frac{Z_i Z_j}{\| r_i - r_j \|^\alpha} & i \neq j\\
            \frac{1}{2} Z_i^{2.4} & i = j
        \end{cases}

    Parameters
    ----------
    numbers : array-like, shape=(n_atoms, )
        The atomic numbers of all the atoms

    coords : array-like, shape=(n_atoms, 3)
        The xyz coordinates of all the atoms (in angstroms)

    alpha : number, default=6
        Some value to exponentiate the distance in the coulomb matrix.

    use_decay : bool, default=False
        This setting defines an extra decay for the values as they get futher
        away from the "central atom". This is to alleviate issues the arise as
        atoms enter or leave the cutoff radius.

    Returns
    -------
    top : array, shape=(n_atoms, n_atoms)
        The coulomb matrix
    """
    top = numpy.outer(numbers, numbers).astype(numpy.float64)
    r = cdist(coords, coords)
    if use_decay:
        other = cdist([coords[0]], coords).reshape(-1)
        r += numpy.add.outer(other, other)

    r **= alpha

    with numpy.errstate(divide='ignore', invalid='ignore'):
        numpy.divide(top, r, top)
    numpy.fill_diagonal(top, 0.5 * numpy.array(numbers) ** 2.4)
    top[top == numpy.Infinity] = 0
    top[numpy.isnan(top)] = 0
    return top def get_element_pairs(elements):
    """
    Extract all the element pairs in a molecule.

    Parameters
    ----------
    elements : list
        All the elements in the molecule

    Returns
    -------
    value : list
        All the element pairs in the molecule
    """
    # This is like computing set(combinations(sorted(elements), 2))
    # We do this because it scales with elements instead of atoms.
    counts = Counter(elements)
    pairs = {}
    order = sorted(counts)
    for i, x in enumerate(order):
        for j, y in enumerate(order):
            if i > j:
                continue
            if x == y and counts[x] < 2:
                continue
            pairs[x, y] = 1
    return list(pairs.keys()) def cosine_decay(R, r_cut=6.):
    r"""
    Compute all the cutoff distances.

    The cutoff is defined as

    .. math::

        f_{R_{c}}(R_{ij}) = \begin{cases}
            0.5 ( \cos( \frac{\pi R_{ij}}{R_c} ) + 1 ), & R_{ij} \le R_c \\
            0,  & otherwise
        \end{cases}


    Parameters
    ----------
    R : array, shape=(N_atoms, N_atoms)
        A distance matrix for all the atoms (scipy.spatial.cdist)

    r_cut : float, default=6.
        The maximum distance allowed for atoms to be considered local to the
        "central atom".

    Returns
    -------
    values : array, shape=(N_atoms, N_atoms)
        The new distance matrix with the cutoff function applied
    """
    values = 0.5 * (numpy.cos(numpy.pi * R / r_cut) + 1)
    values[R > r_cut] = 0
    return values def get_angles(coords):
    r"""
    Get the angles between all triples of coords.

    The resulting values are :math:`[0, \pi]` and all invalid values are NaNs.

    Parameters
    ----------
    coords : numpy.array, shape=(n_atoms, n_dim)
        An array of all the coordinates.

    Returns
    -------
    res : numpy.array, shape=(n_atoms, n_atoms, n_atoms)
        An array the angles of all triples.
    """
    diffs = coords - coords[:, None]
    lengths = numpy.linalg.norm(diffs, axis=2)
    with numpy.errstate(divide='ignore', invalid='ignore'):
        unit = diffs / lengths[:, :, None]
    res = numpy.einsum('ijk,jmk->ijm', unit, unit)
    numpy.clip(res, -1., 1., res)
    numpy.arccos(res, res)
    return res UNKNOWN = 'UNKNOWN'
from builtins import range
from itertools import product
from scipy.spatial.distance import cdist
from .base import BaseFeature, SetMergeMixin, EncodedFeature, FormMixin
from .utils import get_depth_threshold_mask_connections, get_coulomb_matrix
from .utils import get_element_pairs, cosine_decay, get_angles
from .constants import UNKNOWN
import numpy
            local_atoms = ordering[:self.max_occupancy + 1]
            mat = get_coulomb_matrix(numbers[local_atoms],
                                     coords[local_atoms],
                                     alpha=self.alpha,
                                     use_decay=self.use_decay)
            # Take away 1 for the start value
            n = self.max_occupancy - (len(local_atoms) - 1)
            mat = numpy.pad(mat, ((0, n), (0, n)), "constant")
            norm_vals = numpy.linalg.norm(mat, axis=0)
            norm_vals[0] = numpy.inf
            sorting = numpy.argsort(norm_vals)[::-1]
            if self.use_reduced:
                # skip the first value in the diag because it is already in
                # the first row
                diag = numpy.diag(mat)[1:].tolist()
                vectors.append(mat[sorting[0]].tolist() + diag)
            else:
                vectors.append(mat[sorting].flatten())
        return numpy.array(vectors)

    def get_local_coulomb_labels(self):
        idxs = []

        vals = list(range(self.max_occupancy + 1))
        if self.use_reduced:
            for i in vals:
                idxs.append((0, i))
            for i in vals[1:]:
                idxs.append((i, i))
        else:
            for pair in product(vals, vals):
                idxs.append(pair)
        return ['local-coul_%d-%d' % (i, j) for i, j in idxs]


class BehlerParrinello(SetMergeMixin, BaseFeature):
    """
    An implementation of the descriptors used in Behler-Parrinello Neural
    Networks.

    Parameters
    ----------
    input_type : string, default='list'
        Specifies the format the input values will be (must be one of 'list'
        or 'filename').

    n_jobs : int, default=1
        Specifies the number of processes to create when generating the
        features. Positive numbers specify a specifc amount, and numbers less
        than 1 will use the number of cores the computer has.

    r_cut : float, default=6.
        The maximum distance allowed for atoms to be considered local to the
        "central atom".

    r_s : float, default=1.0
        An offset parameter for computing gaussian values between pairwise
        distances.

    eta : float, default=1.0
        A decay parameter for the gaussian distances.

    lambda_ : float, default=1.0
        This value sets the orientation of the cosine function for the angles.
        It should only take values in {-1., 1.}.

    zeta : float, default=1.0
        A decay parameter for the angular terms.

    Attributes
    ----------
    _elements : tuple
        A set of all the elements in the molecules.

    _element_pairs : tuple
        A set of all the element pairs in the molecules.

    References
    ----------
    Behler, J; Parrinello, M. Generalized Neural-Network Representation of
    High-Dimensional Potential-Energy Surfaces. Phys. Rev. Lett. 98, 146401.
    """
    ATTRIBUTES = ("_elements", "_element_pairs")
    LABELS = ("_elements", ("get_chain_labels", "_element_pairs"))

    def __init__(self, input_type='list', n_jobs=1, r_cut=6.0, r_s=1., eta=1.,
                 lambda_=1., zeta=1.):
        super(BehlerParrinello, self).__init__(input_type=input_type,
                                               n_jobs=n_jobs)
        self.r_cut = r_cut
        self.r_s = r_s
        self.eta = eta
        self.lambda_ = lambda_
        self.zeta = zeta
        self._elements = None
        self._element_pairs = None

    def _para_fit(self, X):
        """
        A single instance of the fit procedure.

        This is formulated in a way that the fits can be done completely
        parallel in a map/reduce fashion.

        Parameters
        ----------
        X : object
            An object to use for the fit

        Returns
        -------
        unique_elements : set
            A set of all the elements in the molecule

        pairs : list
            A list of all the element pairs in the molecule
        """
        data = self.convert_input(X)
        # This is just a cheap way to approximate the actual value
        unique_elements = set(data.elements)

        pairs = get_element_pairs(data.elements)