def compute_angles(traj, angle_indices, opt=True):
    """Compute the bond angles between the supplied triplets of indices in each frame of a trajectory.

    Parameters
    ----------
    traj : Trajectory
        An mtraj trajectory.
    angle_indices : np.ndarray, shape=(num_pairs, 2), dtype=int
       Each row gives the indices of three atoms which together make an angle.
    opt : bool, default=True
        Use an optimized native library to calculate distances. Our optimized
        SSE angle calculation implementation is 10-20x faster than the
        (itself optimized) numpy implementation.

    Returns
    -------
    angles : np.ndarray, shape=[n_frames, n_angles], dtype=float
        The angles are in radians
    """
    xyz = ensure_type(traj.xyz, dtype=np.float32, ndim=3, name='traj.xyz', shape=(None, None, 3), warn_on_cast=False)
    triplets = ensure_type(np.asarray(angle_indices), dtype=np.int32, ndim=2, name='angle_indices', shape=(None, 3), warn_on_cast=False)
    if not np.all(np.logical_and(triplets < traj.n_atoms, triplets >= 0)):
        raise ValueError('angle_indices must be between 0 and %d' % traj.n_atoms)

    out = np.zeros((xyz.shape[0], triplets.shape[0]), dtype=np.float32)
    if opt:
        _geometry._angle(xyz, triplets, out)
    else:
        _angle(xyz, triplets, out)
    return out

def compute_dihedrals(trajectory, indices, opt=True):
    """Compute the dihedral angles between the supplied quartets of atoms in each frame in a trajectory.

    Parameters
    ----------
    trajectory : Trajectory
        An mtraj trajectory.
    indices : np.ndarray, shape=(n_dihedrals, 4), dtype=int
        Each row gives the indices of four atoms which together make a
        dihedral angle. The angle is between the planes spanned by the first
        three atoms and the last three atoms, a torsion around the bond
        between the middle two atoms.
    opt : bool, default=True
        Use an optimized native library to calculate angles.

    Returns
    -------
    dihedrals : np.ndarray, shape=(n_frames, n_dihedrals), dtype=float
        The output array gives, in each frame from the trajectory, each of the
        `n_dihedrals` torsion angles. The angles are measured in **radians**.

    """
    xyz = ensure_type(trajectory.xyz, dtype=np.float32, ndim=3, name='traj.xyz', shape=(None, None, 3), warn_on_cast=False)
    quartets = ensure_type(np.asarray(indices), dtype=np.int32, ndim=2, name='indices', shape=(None, 4), warn_on_cast=False)
    if not np.all(np.logical_and(quartets < trajectory.n_atoms, quartets >= 0)):
        raise ValueError('indices must be between 0 and %d' % trajectory.n_atoms)

    out = np.zeros((xyz.shape[0], quartets.shape[0]), dtype=np.float32)
    if opt:
        _geometry._dihedral(xyz, quartets, out)
    else:
        _dihedral(xyz, quartets, out)
    return out

    def write(self, coordinates, topology, time=None, unitcell_vectors=None,
              precision=3):
        """Write one or more frames of a molecular dynamics trajectory to disk
        in the GROMACS GRO format.

        Parameters
        ----------
        coordinates : np.ndarray, dtype=np.float32, shape=(n_frames, n_atoms, 3)
            The cartesian coordinates of each atom, in units of nanometers.
        topology : mdtraj.Topology
            The Topology defining the model to write.
        time : np.ndarray, dtype=float32, shape=(n_frames), optional
            The simulation time corresponding to each frame, in picoseconds.
            If not supplied, the numbers 0..n_frames will be written.
        unitcell_vectors : np.ndarray, dtype=float32, shape=(n_frames, 3, 3), optional
            The periodic box vectors of the simulation in each frame, in nanometers.
        precision : int, optional
            The number of decimal places to print for coordinates. Default is 3
        """
        if not self._open:
            raise ValueError('I/O operation on closed file')
        if not self._mode == 'w':
            raise ValueError('file not opened for writing')

        coordinates = ensure_type(coordinates, dtype=np.float32, ndim=3, name='coordinates', can_be_none=False, warn_on_cast=False)
        time = ensure_type(time, dtype=float, ndim=1, name='time', can_be_none=True, shape=(len(coordinates),), warn_on_cast=False)
        unitcell_vectors = ensure_type(unitcell_vectors, dtype=float, ndim=3, name='unitcell_vectors',
            can_be_none=True, shape=(len(coordinates), 3, 3), warn_on_cast=False)

        for i in range(coordinates.shape[0]):
            frame_time = None if time is None else time[i]
            frame_box = None if unitcell_vectors is None else unitcell_vectors[i]
            self._write_frame(coordinates[i], topology, frame_time, frame_box, precision)

def permute_energies(X, s):
    """Re-order an observable X so that u[i, j, k] correponds to frame i, sampled from state j, evaluated in state k.

    Parameters
    ----------

    X : np.ndarray, shape=(n_iter, n_replicas, n_replicas)
        The observable to permute
    s : np.ndarray, shape=(n_iter, n_replicas), dtype='int'
        The thermodynamic state indices of each replica slot.  s[i, k] is the 
        thermodynamic state index of frame i, replica k.  
    """

    X = ensure_type(X, 'float32', 3, "X")
    n_iter, n_replicas, n_replicas = X.shape
    s = ensure_type(s, "int", 2, "s", shape=(n_iter, n_replicas))
    
    u = np.zeros((n_iter, n_replicas, n_replicas))
    for i, si in enumerate(s):
        mapping = dict(zip(range(n_replicas), si))
        inv_map = {v:k for k, v in mapping.items()}
        si_inv = [inv_map[k] for k in range(n_replicas)]
        u[i] = X[i, si_inv]
    
    return u

def find_closest_contact(traj, group1, group2, frame=0, periodic=True):
    """Find the closest contact between two groups of atoms.

    Given a frame of a Trajectory and two groups of atoms, identify the pair of
    atoms (one from each group) that form the closest contact between the two groups.

    Parameters
    ----------
    traj : Trajectory
        An mtraj trajectory.
    group1 : np.ndarray, shape=(num_atoms), dtype=int
        The indices of atoms in the first group.
    group2 : np.ndarray, shape=(num_atoms), dtype=int
        The indices of atoms in the second group.
    frame : int, default=0
        The frame of the Trajectory to take positions from
    periodic : bool, default=True
        If `periodic` is True and the trajectory contains unitcell
        information, we will compute distances under the minimum image
        convention.

    Returns
    -------
    result : tuple (int, int, float)
         The indices of the two atoms forming the closest contact, and the distance between them.
    """
    xyz = ensure_type(traj.xyz, dtype=np.float32, ndim=3, name='traj.xyz', shape=(None, None, 3), warn_on_cast=False)[frame]
    atoms1 = ensure_type(group1, dtype=np.int32, ndim=1, name='group1', warn_on_cast=False)
    atoms2 = ensure_type(group2, dtype=np.int32, ndim=1, name='group2', warn_on_cast=False)
    if periodic and traj._have_unitcell:
        box = ensure_type(traj.unitcell_vectors, dtype=np.float32, ndim=3, name='unitcell_vectors', shape=(len(traj.xyz), 3, 3),
                          warn_on_cast=False)[frame]
    else:
        box = None
    return _geometry._find_closest_contact(xyz, atoms1, atoms2, box)

    def write(self, xyz, cell_lengths=None):
        """Write one or more frames of data to a mdcrd file

        Parameters
        ----------
        xyz : np.ndarray, shape=(n_frames, n_atoms, 3)
            The cartesian coordinates of the atoms to write. By convention, the
            lengths should be in units of angstroms.
        cell_lengths : np.ndarray, shape=(n_frames, 3), dtype=float32, optional
            The length of the periodic box in each frame, in each direction,
            `a`, `b`, `c`. By convention the lengths should be in units
            of angstroms.
        """
        if not self._mode == 'w':
            raise ValueError('write() is only available when file is opened '
                             'in mode="w"')

        xyz = ensure_type(xyz, np.float32, 3, 'xyz', can_be_none=False,
                shape=(None, None, 3), warn_on_cast=False,
                add_newaxis_on_deficient_ndim=True)
        cell_lengths = ensure_type(cell_lengths, np.float32, 2, 'cell_lengths',
                can_be_none=True, shape=(len(xyz), 3), warn_on_cast=False,
                add_newaxis_on_deficient_ndim=True)

        if self._w_has_box is None:
            # this is the first write()
            self._n_atoms = xyz.shape[1]
            self._fh.write('TITLE : Created by MDTraj with %d atoms\n' % self._n_atoms)

            if cell_lengths is None:
                self._w_has_box = False
            else:
                self._w_has_box = True
        elif self._w_has_box is True:
            if cell_lengths is None:
                raise ValueError('This mdcrd file must contain unitcell '
                                 'information')
        elif self._w_has_box is False:
            if cell_lengths is not None:
                raise ValueError('This mdcrd file must not contain unitcell '
                                 'information')
        else:
            raise RuntimeError()

        for i in range(xyz.shape[0]):
            for j, coord in enumerate(xyz[i].reshape(-1)):
                lfdone = False
                out = "%8.3f" % coord
                if len(out) > 8:
                    raise ValueError('Overflow error')
                self._fh.write(out)
                if (j+1) % 10 == 0:
                    self._fh.write("\n")
                    lfdone = True

            if not lfdone:
                self._fh.write("\n")

            if cell_lengths is not None:
                self._fh.write("%8.3f %8.3f %8.3f\n" % tuple(cell_lengths[i]))

def compute_dihedrals(traj, indices, periodic=True, opt=True):
    """Compute the dihedral angles between the supplied quartets of atoms in each frame in a trajectory.

    Parameters
    ----------
    traj : Trajectory
        An mtraj trajectory.
    indices : np.ndarray, shape=(n_dihedrals, 4), dtype=int
        Each row gives the indices of four atoms which together make a
        dihedral angle. The angle is between the planes spanned by the first
        three atoms and the last three atoms, a torsion around the bond
        between the middle two atoms.
    periodic : bool, default=True
        If `periodic` is True and the trajectory contains unitcell
        information, we will treat dihedrals that cross periodic images
        using the minimum image convention.
    opt : bool, default=True
        Use an optimized native library to calculate angles.

    Returns
    -------
    dihedrals : np.ndarray, shape=(n_frames, n_dihedrals), dtype=float
        The output array gives, in each frame from the trajectory, each of the
        `n_dihedrals` torsion angles. The angles are measured in **radians**.

    """
    xyz = ensure_type(traj.xyz, dtype=np.float32, ndim=3, name='traj.xyz', shape=(None, None, 3), warn_on_cast=False)
    quartets = ensure_type(indices, dtype=np.int32, ndim=2, name='indices', shape=(None, 4), warn_on_cast=False)
    if not np.all(np.logical_and(quartets < traj.n_atoms, quartets >= 0)):
        raise ValueError('indices must be between 0 and %d' % traj.n_atoms)

    if len(quartets) == 0:
        return np.zeros((len(xyz), 0), dtype=np.float32)

    if periodic and traj._have_unitcell:
        if opt and not np.allclose(traj.unitcell_angles, 90):
            warnings.warn('Optimized dihedral calculation does not work for non-orthorhombic '
                          'unit cells and periodic boundary conditions. Falling back to much '
                          'slower pure-Python implementation. Set periodic=False or opt=False '
                          'to disable this message.')
            opt = False

    out = np.zeros((xyz.shape[0], quartets.shape[0]), dtype=np.float32)
    if periodic and traj._have_unitcell:
        box = ensure_type(traj.unitcell_vectors, dtype=np.float32, ndim=3, name='unitcell_vectors', shape=(len(xyz), 3, 3))
        if opt:
            _geometry._dihedral_mic(xyz, quartets, box, out)
            return out
        else:
            _dihedral(traj, quartets, periodic, out)
            return out

    if opt:
        _geometry._dihedral(xyz, quartets, out)
    else:
        _dihedral(traj, quartets, periodic, out)
    return out

def compute_angles(traj, angle_indices, periodic=True, opt=True):
    """Compute the bond angles between the supplied triplets of indices in each frame of a trajectory.

    Parameters
    ----------
    traj : Trajectory
        An mdtraj trajectory.
    angle_indices : np.ndarray, shape=(num_angles, 3), dtype=int
       Each row gives the indices of three atoms which together make an angle.
    periodic : bool, default=True
        If `periodic` is True and the trajectory contains unitcell
        information, we will treat angles that cross periodic images using
        the minimum image convention.
    opt : bool, default=True
        Use an optimized native library to calculate distances. Our optimized
        SSE angle calculation implementation is 10-20x faster than the
        (itself optimized) numpy implementation.

    Returns
    -------
    angles : np.ndarray, shape=[n_frames, n_angles], dtype=float
        The angles are in radians
    """
    xyz = ensure_type(traj.xyz, dtype=np.float32, ndim=3, name='traj.xyz', shape=(None, None, 3), warn_on_cast=False)
    triplets = ensure_type(angle_indices, dtype=np.int32, ndim=2, name='angle_indices', shape=(None, 3), warn_on_cast=False)
    if not np.all(np.logical_and(triplets < traj.n_atoms, triplets >= 0)):
        raise ValueError('angle_indices must be between 0 and %d' % traj.n_atoms)

    if len(triplets) == 0:
        return np.zeros((len(xyz), 0), dtype=np.float32)

    if periodic and traj._have_unitcell:
        if opt and not np.allclose(traj.unitcell_angles, 90):
            warnings.warn('Optimized angle calculation does not work for non-orthorhombic '
                          'unit cells and periodic boundary conditions. Falling back to much '
                          'slower pure-Python implementation. Set periodic=False or opt=False '
                          'to disable this message.')
            opt = False

    out = np.zeros((xyz.shape[0], triplets.shape[0]), dtype=np.float32)
    if periodic is True and traj._have_unitcell:
        box = ensure_type(traj.unitcell_vectors, dtype=np.float32, ndim=3, name='unitcell_vectors', shape=(len(xyz), 3, 3))
        if opt:
            _geometry._angle_mic(xyz, triplets, box, out)
            return out
        else:
            _angle(traj, triplets, periodic, out)
            return out

    if opt:
        _geometry._angle(xyz, triplets, out)
    else:
        _angle(traj, triplets, periodic, out)
    return out

def compute_angles(traj, angle_indices, periodic=True, opt=True):
    """Compute the bond angles between the supplied triplets of indices in each frame of a trajectory.

    Parameters
    ----------
    traj : Trajectory
        An mdtraj trajectory.
    angle_indices : np.ndarray, shape=(num_angles, 3), dtype=int
       Each row gives the indices of three atoms which together make an angle.
    periodic : bool, default=True
        If `periodic` is True and the trajectory contains unitcell
        information, we will treat angles that cross periodic images using
        the minimum image convention.
    opt : bool, default=True
        Use an optimized native library to calculate distances. Our optimized
        SSE angle calculation implementation is 10-20x faster than the
        (itself optimized) numpy implementation.

    Returns
    -------
    angles : np.ndarray, shape=[n_frames, n_angles], dtype=float
        The angles are in radians
    """
    xyz = ensure_type(traj.xyz, dtype=np.float32, ndim=3, name="traj.xyz", shape=(None, None, 3), warn_on_cast=False)
    triplets = ensure_type(
        angle_indices, dtype=np.int32, ndim=2, name="angle_indices", shape=(None, 3), warn_on_cast=False
    )
    if not np.all(np.logical_and(triplets < traj.n_atoms, triplets >= 0)):
        raise ValueError("angle_indices must be between 0 and %d" % traj.n_atoms)

    if len(triplets) == 0:
        return np.zeros((len(xyz), 0), dtype=np.float32)

    out = np.zeros((xyz.shape[0], triplets.shape[0]), dtype=np.float32)
    if periodic is True and traj._have_unitcell:
        box = ensure_type(
            traj.unitcell_vectors, dtype=np.float32, ndim=3, name="unitcell_vectors", shape=(len(xyz), 3, 3)
        )
        if opt:
            orthogonal = np.allclose(traj.unitcell_angles, 90)
            _geometry._angle_mic(xyz, triplets, box.transpose(0, 2, 1).copy(), out, orthogonal)
            return out
        else:
            _angle(traj, triplets, periodic, out)
            return out

    if opt:
        _geometry._angle(xyz, triplets, out)
    else:
        _angle(traj, triplets, periodic, out)
    return out

def validate_input_arrays(predictions, measurements, uncertainties, prior_pops=None):
    """Check input data for correct shape and dtype

    Parameters
    ----------
    predictions : ndarray, shape = (num_frames, num_measurements)
        predictions[j, i] gives the ith observabled predicted at frame j
    measurements : ndarray, shape = (num_measurements)
        measurements[i] gives the ith experimental measurement
    uncertainties : ndarray, shape = (num_measurements)
        uncertainties[i] gives the uncertainty of the ith experiment
    prior_pops : ndarray, shape = (num_frames), optional
        Prior populations of each conformation.  If None, skip.
    
    Notes
    -----
    All inputs must have float64 type and compatible shapes.
    """
    num_frames, num_measurements = predictions.shape

    ensure_type(predictions, np.float64, 2, "predictions")
    ensure_type(measurements, np.float64, 1, "measurements", shape=(num_measurements,))
    ensure_type(uncertainties, np.float64, 1, "uncertainties", shape=(num_measurements,))

    if prior_pops is not None:
        ensure_type(prior_pops, np.float64, 1, "prior_pops", shape=(num_frames,))

def compute_distances(traj, atom_pairs, periodic=True, opt=True):
    """Compute the distances between pairs of atoms in each frame.

    Parameters
    ----------
    traj : Trajectory
        An mtraj trajectory.
    atom_pairs : np.ndarray, shape=(num_pairs, 2), dtype=int
        Each row gives the indices of two atoms involved in the interaction.
    periodic : bool, default=True
        If `periodic` is True and the trajectory contains unitcell
        information, we will compute distances under the minimum image
        convention.
    opt : bool, default=True
        Use an optimized native library to calculate distances. Our optimized
        SSE minimum image convention calculation implementation is over 1000x
        faster than the naive numpy implementation.

    Returns
    -------
    distances : np.ndarray, shape=(n_frames, num_pairs), dtype=float
        The distance, in each frame, between each pair of atoms.
    """
    xyz = ensure_type(traj.xyz, dtype=np.float32, ndim=3, name='traj.xyz', shape=(None, None, 3), warn_on_cast=False)
    pairs = ensure_type(atom_pairs, dtype=np.int32, ndim=2, name='atom_pairs', shape=(None, 2), warn_on_cast=False)
    if not np.all(np.logical_and(pairs < traj.n_atoms, pairs >= 0)):
        raise ValueError('atom_pairs must be between 0 and %d' % traj.n_atoms)

    if len(pairs) == 0:
        return np.zeros((len(xyz), 0), dtype=np.float32)

    if periodic and traj._have_unitcell:
        box = ensure_type(traj.unitcell_vectors, dtype=np.float32, ndim=3, name='unitcell_vectors', shape=(len(xyz), 3, 3),
                          warn_on_cast=False)
        orthogonal = np.allclose(traj.unitcell_angles, 90)
        if opt:
            out = np.empty((xyz.shape[0], pairs.shape[0]), dtype=np.float32)
            _geometry._dist_mic(xyz, pairs, box.transpose(0, 2, 1).copy(), out, orthogonal)
            return out
        else:
            return _distance_mic(xyz, pairs, box.transpose(0, 2, 1), orthogonal)

    # either there are no unitcell vectors or they dont want to use them
    if opt:
        out = np.empty((xyz.shape[0], pairs.shape[0]), dtype=np.float32)
        _geometry._dist(xyz, pairs, out)
        return out
    else:
        return _distance(xyz, pairs)

 def score(self, data):
     """Log-likelihood of sequences under the model
     """
     sequences = [ensure_type(s, dtype=np.float32, ndim=2, name="s") for s in data]
     self.inferrer._sequences = data
     logprob, _ = self.inferrer.do_mslds_estep()
     return logprob

    def write(self, xyz, types=None):
        """Write one or more frames of data to a xyz file.

        Parameters
        ----------
        xyz : np.ndarray, shape=(n_frames, n_atoms, 3)
            The cartesian coordinates of the atoms to write.
        types : np.ndarray, shape(3, )
            The type of each particle.
        """

        if not self._mode == 'w':
            raise ValueError('write() is only available when file is opened '
                             'in mode="w"')

        if not types:
            # Make all particles the same type.
            types = ['X' for _ in xrange(xyz.shape[1])]
        xyz = ensure_type(xyz, np.float32, 3, 'xyz', can_be_none=False,
                        shape=(None, None, 3), warn_on_cast=False,
                        add_newaxis_on_deficient_ndim=True)
        in_units_of(xyz, 'nanometers', self.distance_unit, inplace=True)

        for i in range(xyz.shape[0]):
            self._fh.write('{0}\n'.format(xyz.shape[1]))
            self._fh.write("Created with MDTraj {0}, {1}\n".format(version, str(date.today())))

            for j, coord in enumerate(xyz[i]):
                self._fh.write('{0} {1:8.3f} {2:8.3f} {3:8.3f}\n'.format(
                    types[j], coord[0], coord[1], coord[2]))

def shrake_rupley(traj, probe_radius=0.14, n_sphere_points=960):
    """Compute the solvent accessible surface area of each atom in each simulation frame.

    Parameters
    ----------
    traj : Trajectory
        An mtraj trajectory.
    probe_radius : float, optional
        The radius of the probe, in nm.
    n_sphere_pts : int, optional
        The number of points representing the surface of each atom, higher
        values leads to more accuracy.

    Returns
    -------
    areas : np.array, shape=(n_frames, n_atoms)
        The accessible surface area of each atom in every frame

    Notes
    -----
    This code implements the Shrake and Rupley algorithm, with the Golden
    Section Spiral algorithm to generate the sphere points. The basic idea
    is to great a mesh of points representing the surface of each atom
    (at a distance of the van der waals radius plus the probe
    radius from the nuclei), and then count the number of such mesh points
    that are on the molecular surface -- i.e. not within the radius of another
    atom. Assuming that the points are evenly distributed, the number of points
    is directly proportional to the accessible surface area (its just 4*pi*r^2
    time the fraction of the points that are accessible).

    There are a number of different ways to generate the points on the sphere --
    possibly the best way would be to do a little "molecular dyanmics" : put the
    points on the sphere, and then run MD where all the points repel one another
    and wait for them to get to an energy minimum. But that sounds expensive.

    This code uses the golden section spiral algorithm
    (picture at http://xsisupport.com/2012/02/25/evenly-distributing-points-on-a-sphere-with-the-golden-sectionspiral/)
    where you make this spiral that traces out the unit sphere and then put points
    down equidistant along the spiral. It's cheap, but not perfect.

    The gromacs utility g_sas uses a slightly different algorithm for generating
    points on the sphere, which is based on an icosahedral tesselation.
    roughly, the icosahedral tesselation works something like this
    http://www.ziyan.info/2008/11/sphere-tessellation-using-icosahedron.html

    References
    ----------
    .. [1] Shrake, A; Rupley, JA. (1973) J Mol Biol 79 (2): 351--71.
    """
    if not _geometry._processor_supports_sse41():
        raise RuntimeError('This CPU does not support the required instruction set (SSE4.1)')

    xyz = ensure_type(traj.xyz, dtype=np.float32, ndim=3, name='traj.xyz', shape=(None, None, 3), warn_on_cast=False)
    out = np.zeros((xyz.shape[0], xyz.shape[1]), dtype=np.float32)
    atom_radii = [_ATOMIC_RADII[atom.element.symbol] for atom in traj.topology.atoms]
    radii = np.array(atom_radii, np.float32) + probe_radius

    _geometry._sasa(xyz, radii, int(n_sphere_points), out)

    return out

    def _init(self, sequences, init_params):
        """Find initial means(hot start)"""
        sequences = [ensure_type(s, dtype=np.float32, ndim=2, name='s', warn_on_cast=False)
                     for s in sequences]
        self._impl._sequences = sequences

        if self.n_hotstart == 'all':
            small_dataset = np.vstack(sequences)
        else:
            small_dataset = np.vstack(sequences[0:min(len(sequences), self.n_hotstart)])

        if self.init_algo == "GMM" and ("m" in init_params or "v" in init_params):
            mixture = sklearn.mixture.GMM(self.n_states, n_init=1, random_state=self.random_state)
            mixture.fit(small_dataset)
            if "m" in init_params:
                self.means_ = mixture.means_
            if "v" in init_params:
                self.vars_ = mixture.covars_
        else:
            if 'm' in init_params:
                with warnings.catch_warnings():
                    warnings.simplefilter("ignore")
                    self.means_ = cluster.KMeans(
                        n_clusters=self.n_states, n_init=1, init='random',
                        n_jobs=self.n_jobs, random_state=self.random_state).fit(
                        small_dataset).cluster_centers_
            if 'v' in init_params:
                self.vars_ = np.vstack([np.var(small_dataset, axis=0)] * self.n_states)
        if 't' in init_params:
            transmat_ = np.empty((self.n_states, self.n_states))
            transmat_.fill(1.0 / self.n_states)
            self.transmat_ = transmat_
            self.populations_ = np.ones(self.n_states) / self.n_states

    def select_pairs(self, selection1=None, selection2=None):
        """Generate unique pairs of atom indices.

        If a selecton is a string, it will be resolved using the atom selection
        DSL, otherwise it is expected to be an array of atom indices.

        Parameters
        ----------
        selection1 : str or array-like, shape=(n_indices, ), dtype=int
            A selection for `select()` or an array of atom indices.
        selection2 : str or array-like, shape=(n_indices, ), dtype=int
            A selection for `select()` or an array of atom indices.

        Returns
        -------
        pairs : array-like, shape=(n_pairs, 2), dtype=int
            Each row gives the indices of two atoms.

        """
        # Resolve selections using the atom selection DSL...
        if isinstance(selection1, string_types):
            a_indices = self.select(selection1)
        else:  # ...or use a provided array of indices.
            a_indices = ensure_type(selection1, dtype=np.int32, ndim=1,
                                    name='a_indices', warn_on_cast=False)
        if isinstance(selection2, string_types):
            b_indices = self.select(selection2)
        else:
            b_indices = ensure_type(selection2, dtype=np.int32, ndim=1,
                                    name='b_indices', warn_on_cast=False)
        a_indices.sort()
        b_indices.sort()

        # Create unique pairs from the indices.
        if np.array_equal(a_indices, b_indices):
            # This is more efficient and memory friendly by removing the
            # intermediate set creation required in the case below.
            pairs = np.fromiter(itertools.chain.from_iterable(itertools.combinations(a_indices, 2)),
                                dtype=np.int32, count=len(a_indices) * (len(a_indices) - 1))
            pairs = np.vstack((pairs[::2], pairs[1::2])).T
        else:
            pairs = np.array(list(set(
                (a, b) if a > b else (b, a)
                for a, b in itertools.product(a_indices, b_indices)
                if a != b)),
                             dtype=np.int32)
        return pairs

    def select_pairs(self, selection1=None, selection2=None):
        """Generate unique pairs of atom indices.

        If a selecton is a string, it will be resolved using the atom selection
        DSL, otherwise it is expected to be an array of atom indices.

        Parameters
        ----------
        selection1 : str or array-like, shape=(n_indices, ), dtype=int
            A selection for `select()` or an array of atom indices.
        selection2 : str or array-like, shape=(n_indices, ), dtype=int
            A selection for `select()` or an array of atom indices.

        Returns
        -------
        pairs : array-like, shape=(n_pairs, 2), dtype=int
            Each row gives the indices of two atoms.

        """
        # Resolve selections using the atom selection DSL...
        if isinstance(selection1, string_types):
            a_indices = self.select(selection1)
        else:  # ...or use a provided array of indices.
            a_indices = ensure_type(selection1, dtype=np.int32, ndim=1,
                                    name='a_indices', warn_on_cast=False)
        if isinstance(selection2, string_types):
            b_indices = self.select(selection2)
        else:
            b_indices = ensure_type(selection2, dtype=np.int32, ndim=1,
                                    name='b_indices', warn_on_cast=False)
        a_indices.sort()
        b_indices.sort()

        # Create unique pairs from the indices.
        # In the cases where a_indices and b_indices are identical or mutually
        # exclusive, we can utilize a more efficient and memory friendly
        # approach by removing the intermediate set creation required in
        # the general case.
        if np.array_equal(a_indices, b_indices):
            pairs = self._unique_pairs_equal(a_indices)
        elif len(np.intersect1d(a_indices, b_indices)) == 0:
            pairs = self._unique_pairs_mutually_exclusive(a_indices, b_indices)
        else:
            pairs = self._unique_pairs(a_indices, b_indices)
        return pairs

def test_ensure_type_25():
    with warnings.catch_warnings(record=True) as w:
        warnings.simplefilter("always")

        val = ensure_type(a, np.float64, 1, '', length=10, warn_on_cast=False)

        assert val.dtype == np.float64
        assert a.dtype == np.float32  # a should not be changed
        assert len(w) == 0  # no warning since we set warn_on_cast to False

def test_ensure_type_2():
    with warnings.catch_warnings(record=True) as w:
        warnings.simplefilter("always")

        val = ensure_type(a, np.float64, 1, '', length=10)

        assert val.dtype == np.float64
        assert a.dtype == np.float32  # a should not be changed
        assert len(w) == 1
        assert issubclass(w[-1].category, TypeCastPerformanceWarning)

def compute_displacements(traj, atom_pairs, periodic=True, opt=True):
    """Compute the displacement vector between pairs of atoms in each frame of a trajectory.

    Parameters
    ----------
    traj : Trajectory
        Trajectory to compute distances in
    atom_pairs : np.ndarray, shape[num_pairs, 2], dtype=int
        Each row gives the indices of two atoms.
    periodic : bool, default=True
        If `periodic` is True and the trajectory contains unitcell
        information, we will compute distances under the minimum image
        convention.
    opt : bool, default=True
        Use an optimized native library to calculate distances. Our
        optimized minimum image convention calculation implementation is
        over 1000x faster than the naive numpy implementation.

    Returns
    -------
    displacements : np.ndarray, shape=[n_frames, n_pairs, 3], dtype=float32
         The displacememt vector, in each frame, between each pair of atoms.
    """
    xyz = ensure_type(traj.xyz, dtype=np.float32, ndim=3, name='traj.xyz', shape=(None, None, 3))
    pairs = ensure_type(np.asarray(atom_pairs), dtype=np.int32, ndim=2, name='atom_pairs', shape=(None, 2))
    if not np.all(np.logical_and(pairs < traj.n_atoms, pairs >= 0)):
        raise ValueError('atom_pairs must be between 0 and %d' % traj.n_atoms)

    if periodic is True and traj._have_unitcell:
        box = ensure_type(traj.unitcell_vectors, dtype=np.float32, ndim=3, name='unitcell_vectors', shape=(len(xyz), 3, 3))
        if opt and _geometry._processor_supports_sse41():
            out = np.empty((xyz.shape[0], pairs.shape[0], 3), dtype=np.float32)
            _geometry._dist_mic_displacement(xyz, pairs, box, out)
            return out
        else:
            return _displacement_mic(xyz, pairs, box)

    # either there are no unitcell vectors or they dont want to use them
    if opt and _geometry._processor_supports_sse41():
        out = np.empty((xyz.shape[0], pairs.shape[0], 3), dtype=np.float32)
        _geometry._dist_displacement(xyz, pairs, out)
        return out
    return _displacement(xyz, pairs)