def vRaman(x,omega=1.0,delta=0.0,epsilon=0.048,phi=4.0/3.0):
    x=np.array(x)    
    s=1
    v=v=np.einsum('i,jk->ijk',omega*np.exp(1.0j*2.0*phi*x),Fx(s)*np.sqrt(2.0)/2.0)
    v=np.array([np.triu(v[i])+np.conjugate(np.triu(v[i],1)).transpose() for i in range(x.size)])
    v+=np.array([np.diag([epsilon+delta,0.0,epsilon-delta])]*x.size)
    return v

 def pressure(self):
     ps = 1.
     sig = 1.
     eps = 1.
     d = 2.
     N = np.size(self.c.x)
     
     dx = self.c.dx()
     dy = self.c.dy()
     dz = self.c.dz()
     dr2 = self.c.dr2()
     
     r_mag = (dx ** 2 + dy ** 2 + dz ** 2)
     r_mag = np.nan_to_num(r_mag)
     
     px = dx * (24 * eps / r_mag) * ((2 * (sig**12 / r_mag**6)) - (sig**6 / r_mag**3))
     px = np.nan_to_num(px)
     px = np.triu(px)
     py = dy * (24 * eps / r_mag) * ((2 * (sig**12 / r_mag**6)) - (sig**6 / r_mag**3))
     py = np.nan_to_num(py)
     py = np.triu(py)
     pz = dz * (24 * eps / r_mag) * ((2 * (sig**12 / r_mag**6)) - (sig**6 / r_mag**3))
     pz = np.nan_to_num(pz)
     pz = np.triu(pz)
     pt = px + py + pz
     return 1 / (d * N * self.ke()) * np.sum(pt)

    def test_dynamic_programming_logic(self):
        # Test for the dynamic programming part
        # This test is directly taken from Cormen page 376.
        arrays = [np.random.random((30, 35)),
                  np.random.random((35, 15)),
                  np.random.random((15, 5)),
                  np.random.random((5, 10)),
                  np.random.random((10, 20)),
                  np.random.random((20, 25))]
        m_expected = np.array([[0., 15750., 7875., 9375., 11875., 15125.],
                               [0.,     0., 2625., 4375.,  7125., 10500.],
                               [0.,     0.,    0.,  750.,  2500.,  5375.],
                               [0.,     0.,    0.,    0.,  1000.,  3500.],
                               [0.,     0.,    0.,    0.,     0.,  5000.],
                               [0.,     0.,    0.,    0.,     0.,     0.]])
        s_expected = np.array([[0,  1,  1,  3,  3,  3],
                               [0,  0,  2,  3,  3,  3],
                               [0,  0,  0,  3,  3,  3],
                               [0,  0,  0,  0,  4,  5],
                               [0,  0,  0,  0,  0,  5],
                               [0,  0,  0,  0,  0,  0]], dtype=np.int)
        s_expected -= 1  # Cormen uses 1-based index, python does not.

        s, m = _multi_dot_matrix_chain_order(arrays, return_costs=True)

        # Only the upper triangular part (without the diagonal) is interesting.
        assert_almost_equal(np.triu(s[:-1, 1:]),
                            np.triu(s_expected[:-1, 1:]))
        assert_almost_equal(np.triu(m), np.triu(m_expected))

    def get_representational_similarity(self, corr='spearman'):
        logging.debug("Calculating representational similarity")

        response = self.response[:, 1:, :self.numbercells, 0] # orientation x freq x phase x cell, no blank
        response = response.reshape(self.number_ori * (self.number_tf-1), self.numbercells)
        Nstim, N = response.shape

        rep_sim = np.zeros((Nstim, Nstim))
        rep_sim_p = np.empty((Nstim, Nstim))
        if corr == 'pearson':
            for i in range(Nstim):
                for j in range(i, Nstim): # matrix is symmetric
                    rep_sim[i, j], rep_sim_p[i, j] = st.pearsonr(response[i], response[j])

        elif corr == 'spearman':
            for i in range(Nstim):
                for j in range(i, Nstim): # matrix is symmetric
                    rep_sim[i, j], rep_sim_p[i, j] = st.spearmanr(response[i], response[j])

        else:
            raise Exception('correlation should be pearson or spearman')

        rep_sim = np.triu(rep_sim) + np.triu(rep_sim, 1).T # fill in lower triangle
        rep_sim_p = np.triu(rep_sim_p) + np.triu(rep_sim_p, 1).T  # fill in lower triangle

        return rep_sim, rep_sim_p

    def test_tpsv(self):
        seed(1234)
        for ind, dtype in enumerate(DTYPES):
            n = 10
            x = rand(n).astype(dtype)
            # Upper triangular array
            A = triu(rand(n, n)) if ind < 2 else triu(rand(n, n)+rand(n, n)*1j)
            A += eye(n)
            # Form the packed storage
            c, r = tril_indices(n)
            Ap = A[r, c]
            func, = get_blas_funcs(('tpsv',), dtype=dtype)

            y1 = func(n=n, ap=Ap, x=x)
            y2 = solve(A, x)
            assert_array_almost_equal(y1, y2)

            y1 = func(n=n, ap=Ap, x=x, diag=1)
            A[arange(n), arange(n)] = dtype(1)
            y2 = solve(A, x)
            assert_array_almost_equal(y1, y2)

            y1 = func(n=n, ap=Ap, x=x, diag=1, trans=1)
            y2 = solve(A.T, x)
            assert_array_almost_equal(y1, y2)

            y1 = func(n=n, ap=Ap, x=x, diag=1, trans=2)
            y2 = solve(A.conj().T, x)
            assert_array_almost_equal(y1, y2)

    def __init__(self, endmembers, alphas, energy_interaction, volume_interaction=None, entropy_interaction=None):

        self.n_endmembers = len(endmembers)

        # Create array of van Laar parameters
        self.alphas = np.array(alphas)

        # Create 2D arrays of interaction parameters
        self.We = np.triu(2. / (self.alphas[:, np.newaxis] + self.alphas), 1)
        self.We[np.triu_indices(self.n_endmembers, 1)] *= np.array([i for row in energy_interaction
                                                                for i in row])

        if entropy_interaction is not None:
            self.Ws = np.triu(2. / (self.alphas[:, np.newaxis] + self.alphas), 1)
            self.Ws[np.triu_indices(self.n_endmembers, 1)] *= np.array([i for row in entropy_interaction
                                                                        for i in row])
        else:
            self.Ws = np.zeros((self.n_endmembers, self.n_endmembers))

        if volume_interaction is not None:
            self.Wv = np.triu(2. / (self.alphas[:, np.newaxis] + self.alphas), 1)
            self.Wv[np.triu_indices(self.n_endmembers, 1)] *= np.array([i for row in volume_interaction
                                                                        for i in row])
        else:
            self.Wv = np.zeros((self.n_endmembers, self.n_endmembers))


        # initialize ideal solution model
        IdealSolution.__init__(self, endmembers)

    def compute_distances_and_pairs(self, pdb_file, nr_contacts=None, nr_noncontacts=None):
        #distance and contacts
        self.features['pair']['Cbdist'] = pdb.distance_map(pdb_file, self.L)

        #mask positions that have too many gaps
        gap_freq = 1 - (self.Ni / self.neff)
        highly_gapped_pos = np.where(gap_freq > self.max_gap_percentage)[0]
        self.features['pair']['Cbdist'][:,highly_gapped_pos] = np.nan
        self.features['pair']['Cbdist'][highly_gapped_pos, :] = np.nan

        #if there are unresolved residues, there will be nan in the distance_map
        with np.errstate(invalid='ignore'):
            self.features['pair']['contact'] = (self.features['pair']['Cbdist'] <= self.contact_threshold) * 1
            self.features['pair']['nocontact'] = (self.features['pair']['Cbdist'] > self.non_contact_threshold) * 1

        indices_contact = np.where(np.triu(self.features['pair']['contact'], k=self.seq_separation))
        indices_contact = tuple(shuffle(indices_contact[0],indices_contact[1], random_state=0))
        if nr_contacts:
            indices_contact = indices_contact[0][:nr_contacts], indices_contact[1][:nr_contacts]

        indices_nocontact = np.where(np.triu(self.features['pair']['nocontact'], k=self.seq_separation))
        indices_nocontact = tuple(shuffle(indices_nocontact[0],indices_nocontact[1], random_state=0))
        if nr_noncontacts:
            indices_nocontact = indices_nocontact[0][:nr_noncontacts], indices_nocontact[1][:nr_noncontacts]


        #update indices of i<j for only relevant pairs
        self.ij_ind_upper = np.array(list(indices_contact[0]) + list(indices_nocontact[0])), np.array(list(indices_contact[1]) + list(indices_nocontact[1]))

	def __init__(self, num_visibles, num_hiddens):
		"""
		Initializes the parameters of the SemiRBM.

		@type  num_visibles: integer
		@param num_visibles: number of visible units
		@type  num_hiddens:  integer
		@param num_hiddens:  number of hidden units
		"""

		AbstractBM.__init__(self, num_visibles, num_hiddens)

		# additional hyperparameters
		self.learning_rate_lateral = 0.01
		self.momentum_lateral = 0.5
		self.weight_decay_lateral = 0.

		self.damping = 0.2
		self.num_lateral_updates = 20
		self.sampling_method = AbstractBM.MF

		# additional parameters
		self.L = np.matrix(np.random.randn(num_visibles, num_visibles)) / num_visibles / 200
		self.L = np.triu(self.L) + np.triu(self.L).T - 2. * np.diag(np.diag(self.L))
		self.dL = np.zeros_like(self.L)

def calculate2_pseudoV(pred, truth, rnd=0.01, full_matrix=True, sym=False):
    if full_matrix:
        pred_cp = pred
        truth_cp = truth
    else: # make matrix upper triangular
        pred_cp = np.triu(pred)
        truth_cp = np.triu(truth)

    # Avoid dividing by zero by rounding everything less than rnd up to rnd
    # Note: it is ok to do this after making the matrix upper triangular
    # since the bottom triangle of the matrix will not affect the score

    size = np.array(pred_cp.shape)[1]
    res = 0 # result to be returned

    # do one row at a time to reduce memory usage
    for x in xrange(size):
        # (1 - rnd) will cast the pred_cp/truth_cp matrices automatically if they are int8
        pred_row = (1 - rnd) * pred_cp[x, ] + rnd
        truth_row = (1 - rnd) * truth_cp[x, ] + rnd

        pred_row /= np.sum(pred_row)
        truth_row /= np.sum(truth_row)
        if sym:
            res += np.sum(truth_row * np.log(truth_row/pred_row)) + np.sum(pred_row * np.log(pred_row/truth_row))
        else:
            res += np.sum(truth_row * np.log(truth_row/pred_row))
    return res

def vRaman(x,omega=1.0,delta=0.0,epsilon=0.048,phi=4.0/3.0):
    x=np.array(x)    
    s=1
    v=np.outer(omega*np.exp(1.0j*2.0*phi*x),Fx(s)*np.sqrt(2.0)/2.0).reshape(x.size,2*s+1,2*s+1)
    v=sLA.block_diag(*[np.triu(v[i])+np.conjugate(np.triu(v[i],1)).transpose() for i in range(x.size)])
    v+=sLA.block_diag(*[np.diag([epsilon-delta,0.0,epsilon+delta])]*x.size)
    return v

 def test_compute_modes(self):
     ws = N.identity(self.num_states)
     tol = 1e-6
     weights_full = N.mat(N.random.random((self.num_states, self.num_states)))
     weights_full = N.triu(weights_full) + N.triu(weights_full, 1).H
     weights_full = weights_full*weights_full
     weights_diag = N.random.random(self.num_states)
     weights_list = [None, weights_diag, weights_full]
     vec_array = N.random.random((self.num_states, self.num_vecs))
     for weights in weights_list:
         IP = VectorSpaceMatrices(weights=weights).compute_inner_product_mat
         correlation_mat_true = IP(vec_array, vec_array)
         eigen_vals_true, eigen_vecs_true = util.eigh(correlation_mat_true)
         build_coeff_mat_true = eigen_vecs_true * N.mat(N.diag(
             eigen_vals_true**-0.5))
         modes_true = vec_array.dot(build_coeff_mat_true)
         
         modes, eigen_vals, eigen_vecs, correlation_mat = \
             compute_POD_matrices_snaps_method(vec_array, self.mode_indices, 
             inner_product_weights=weights, return_all=True)
         
         N.testing.assert_allclose(eigen_vals, eigen_vals_true, rtol=tol)
         N.testing.assert_allclose(eigen_vecs, eigen_vecs_true)
         N.testing.assert_allclose(correlation_mat, correlation_mat_true)
         N.testing.assert_allclose(modes, modes_true[:,self.mode_indices])
                     
         modes, eigen_vals, eigen_vecs = \
             compute_POD_matrices_direct_method(vec_array, self.mode_indices, 
             inner_product_weights=weights, return_all=True)
         
         N.testing.assert_allclose(eigen_vals, eigen_vals_true)
         N.testing.assert_allclose(N.abs(eigen_vecs), N.abs(eigen_vecs_true))
         N.testing.assert_allclose(N.abs(modes), N.abs(modes_true[:,self.mode_indices]))

    def prob(self, state):
        """
        Calculates the probability of a given state using the exponential
        family parameterization learned from a prior dataset. Note that
        we are assuming A(theta) = 0 (i.e., the graph is triangulated).
        Also note the similarity to calc_mu.
        """
        n1 = state[0,0:self.num_sites]
        n2 = state[0,self.num_sites:]
        mu_s = state
        s = np.zeros((self.num_nodes, self.num_nodes))
        mu_st11 = np.matrix(s)
        mu_st10 = np.matrix(s)
        mu_st01 = np.matrix(s)
        mu_st00 = np.matrix(s)

        """
        Calculate the edges from n->n'.
        This is the upper-right quadrant.
        """	
        upright = self.calc_mu_quadrant(n1, n2, False)
        mu_st11[0:self.num_sites,self.num_sites:] += upright[3] # nn11
        mu_st10[0:self.num_sites,self.num_sites:] += upright[2] # nn10
        mu_st01[0:self.num_sites,self.num_sites:] += upright[1] # nn01
        mu_st00[0:self.num_sites,self.num_sites:] += upright[0] # nn00
        """
        Calculate the edges from n'->n'.
        This is the lower-right quadrant.
        """
        lowright = self.calc_mu_quadrant(n2, n2, True)
        mu_st11[self.num_sites:,self.num_sites:] += lowright[3] # nn11
        mu_st10[self.num_sites:,self.num_sites:] += lowright[2] # nn10
        mu_st01[self.num_sites:,self.num_sites:] += lowright[1] # nn01
        mu_st00[self.num_sites:,self.num_sites:] += lowright[0] # nn00
        """
        We only want the upper triangle, since otherwise we will double
        count the n'->n' edges.
        """
        mu_st11 = np.triu(mu_st11)
        mu_st10 = np.triu(mu_st10)
        mu_st01 = np.triu(mu_st01)
        mu_st00 = np.triu(mu_st00)
        """
        Sum(x_s)
        """
        sum_s = mu_s * self.nodes[0].T
        sum_s += (mu_s - 1) * -1 * self.nodes[1].T
        sum_s = sum_s[0,0]
        """
        Sum(x_s,x_t)
        """
        sum_st00 = np.sum(np.multiply(mu_st00, self.edges[0]))
        sum_st01 = np.sum(np.multiply(mu_st01, self.edges[1]))
        sum_st10 = np.sum(np.multiply(mu_st10, self.edges[2]))
        sum_st11 = np.sum(np.multiply(mu_st11, self.edges[3]))
        result = sum_s + sum_st00 + sum_st01 + sum_st10 + sum_st11
        #print "S: {0} ST[00]: {1} ST[01]: {2} ST[10]: {3} ST[11]: {4}".format(sum_s, sum_st00, sum_st01, sum_st10, sum_st11)
        #print "Prob: {0}".format(result)
        #pseudo-likelihood
        return result

def make_connections(n,density=0.25):
	"""
	This function will return a random adjacency matrix of size
	n x n. You read the matrix like this:
	
	if matrix[2,7] = 1, then cities '2' and '7' are connected.
	if matrix[2,7] = 0, then the cities are _not_ connected.
	
	:param n: number of cities
	:param density: controls the ratio of 1s to 0s in the matrix
	
	:returns: an n x n adjacency matrix
	"""
	
	import networkx
	
	# Generate a random adjacency matrix and use it to build a networkx graph
	a=numpy.int32(numpy.triu((numpy.random.random_sample(size=(n,n))<density)))
	G=networkx.from_numpy_matrix(a)
	
	# If the network is 'not connected' (i.e., there are isolated nodes)
	# generate a new one. Keep doing this until we get a connected one.
	# Yes, there are more elegant ways to do this, but I'm demonstrating
	# while loops!
	while not networkx.is_connected(G):
		a=numpy.int32(numpy.triu((numpy.random.random_sample(size=(n,n))<density)))
		G=networkx.from_numpy_matrix(a)
	
	# Cities should be connected to themselves.
	numpy.fill_diagonal(a,1)
	
	return a + numpy.triu(a,1).T

    def __init__(self, rng, matches_vec, batch_size,
            sample_diff_every_epoch=True, n_same_pairs=None):
        """
        If `n_same_pairs` is given, this number of same pairs is sampled,
        otherwise all same pairs are used.
        """
        self.rng = rng
        self.matches_vec = matches_vec
        self.batch_size = batch_size

        if n_same_pairs is None:
            # Use all pairs
            I, J = np.where(np.triu(distance.squareform(matches_vec)))  # indices of same pairs
        else:
            # Sample same pairs
            n_pairs = min(n_same_pairs, len(np.where(matches_vec == True)[0]))
            same_sample = self.rng.choice(
                np.where(matches_vec == True)[0], size=n_pairs, replace=False
                )
            same_vec = np.zeros(self.matches_vec.shape[0], dtype=np.bool)
            same_vec[same_sample] = True
            I, J = np.where(np.triu(distance.squareform(same_vec)))

        self.x1_same_indices = []
        self.x2_same_indices = []
        for i, j in zip(I, J):
            self.x1_same_indices.append(i)
            self.x2_same_indices.append(j)

        if not sample_diff_every_epoch:
            self.x1_diff_indices, self.x2_diff_indices = self._sample_diff_pairs()
        self.sample_diff_every_epoch = sample_diff_every_epoch

    def __call__(self, container, sim_wide_params):
        distance_matrices = container.get_distance_matrices()
        # LJ
        accelerations, pe = moldyn.lennardJonesForce(
                distance_matrices,
                radius=self.radius,
                epsilon=self.epsilon,
                cutoff_dist=self.cutoff_dist,
                anchor_ixs = sim_wide_params.anchor_ixs)
        self.pe = pe  # should LJ func calc this differently since we are doing a cutoff?
        # Damping
        dx, dy, dr = distance_matrices  # KLUDGE
        ixs_not_touching = dr > self.cutoff_dist
        self.ixs_not_touching = ixs_not_touching
        dvx, dvy, dvr = moldyn.get_distance_matrices(container.velocities)
##        if np.any(container.velocities > 1):
##            pass
##        if container.time > 1:
##            pass
        dr[dr == 0] = 0.1  # prevent / 0 errors. Will this cause problems?
        accel = self.damping_magnitude * (dvx*dx + dvy*dy)/dr
        accel[ixs_not_touching] = 0
        accelerations[:, 0] += np.sum(np.triu(accel * dx/dr), axis=0)
        accelerations[:, 1] += np.sum(np.triu(accel * dy/dr), axis=0)
        return accelerations

    def pstep(self, level=0):
        """Propagates the momenta for half a time step."""

        dt = self.pdt[level]
        dt2 = dt**2
        dt3 = dt**3 / 3.0
        m = dstrip(self.beads.m3)[0].reshape(self.beads.natoms, 3)

        hh0 = np.dot(self.cell.h, self.h0.ih)
        pi_ext = np.dot(hh0, np.dot(self.stressext, hh0.T)) * self.h0.V / self.cell.V
        L = np.diag([3, 2, 1])

        stress = dstrip(self.stress_mts(level))
        self.p += dt * (self.cell.V * np.triu(stress))

        # integerates the kinetic part of the stress with the force at the inner-most level.
        if(level == self.nmtslevels - 1):

            self.p += dt * (self.cell.V * np.triu(-self.beads.nbeads * pi_ext) + Constants.kb * self.temp * L)

            pc = dstrip(self.beads.pc).reshape(self.beads.natoms, 3)
            fc = np.sum(dstrip(self.forces.forces_mts(level)), axis=0).reshape(self.beads.natoms, 3) / self.beads.nbeads
            fcTonm = (fc / dstrip(self.beads.m3)[0].reshape(self.beads.natoms, 3)).T

            self.p += np.triu(dt2 * np.dot(fcTonm, pc) + dt3 * np.dot(fcTonm, fc)) * self.beads.nbeads

def test_tzrzf():
    """
    This test performs an RZ decomposition in which an m x n upper trapezoidal
    array M (m <= n) is factorized as M = [R 0] * Z where R is upper triangular
    and Z is unitary.
    """
    seed(1234)
    m, n = 10, 15
    for ind, dtype in enumerate(DTYPES):
        tzrzf, tzrzf_lw = get_lapack_funcs(('tzrzf', 'tzrzf_lwork'),
                                           dtype=dtype)
        lwork = _compute_lwork(tzrzf_lw, m, n)

        if ind < 2:
            A = triu(rand(m, n).astype(dtype))
        else:
            A = triu((rand(m, n) + rand(m, n)*1j).astype(dtype))

        # assert wrong shape arg, f2py returns generic error
        assert_raises(Exception, tzrzf, A.T)
        rz, tau, info = tzrzf(A, lwork=lwork)
        # Check success
        assert_(info == 0)

        # Get Z manually for comparison
        R = np.hstack((rz[:, :m], np.zeros((m, n-m), dtype=dtype)))
        V = np.hstack((np.eye(m, dtype=dtype), rz[:, m:]))
        Id = np.eye(n, dtype=dtype)
        ref = [Id-tau[x]*V[[x], :].T.dot(V[[x], :].conj()) for x in range(m)]
        Z = reduce(np.dot, ref)
        assert_allclose(R.dot(Z) - A, zeros_like(A, dtype=dtype),
                        atol=10*np.spacing(dtype(1.0).real), rtol=0.)

 def _getImageDirection(self, threshold=None):
     """
     _getImageDirection(threshold=None)
     
     Returns the direction of the sequence of avalanches as: 
     "Top_to_bottom","Left_to_right", "Bottom_to_top","Right_to_left"
     
     Parameters:
     ----------------
     threshold : int
         Minimum value of the gray level change to conisider
         a pixel as part of an avalanche (i.e. it is switched)
     """
     # Top, left, bottom, rigth
     imageDirections=["Top_to_bottom","Left_to_right", "Bottom_to_top","Right_to_left"]
     # Check if color Image is available
     if not self._isColorImage:
         self._isColorImageDone(ask=False)        
     switchTimesMasked = self._switchTimes2D
     pixelsUnderMasks = []
     # first identify first 10 avalanches of whole image
     firstAvsList = np.unique(self._switchTimes2D)[:11]
     # Prepare the mask
     m = np.ones((self.dimX, self.dimY))
     # Top mask
     mask = np.rot90(np.triu(m)) * np.triu(m)
     top = switchTimesMasked * mask
     pixelsUnderMasks.append(sum([np.sum(top==elem) for elem in firstAvsList]))
     # Now we need to rotate the mask
     for i in range(3):
         mask = np.rot90(mask)
         top = switchTimesMasked * mask
         pixelsUnderMasks.append(sum([np.sum(top==elem) for elem in firstAvsList]))
     max_in_mask = scipy.array(pixelsUnderMasks).argmax()
     return imageDirections[max_in_mask]

  def test_iris(self):

    # Generate full set of constraints for comparison with reference implementation
    mask = (self.iris_labels[None] == self.iris_labels[:, None])
    a, b = np.nonzero(np.triu(mask, k=1))
    c, d = np.nonzero(np.triu(~mask, k=1))

    # Full metric
    mmc = MMC(convergence_threshold=0.01)
    mmc.fit(self.iris_points, [a, b, c, d])
    expected = [[+0.00046504, +0.00083371, -0.00111959, -0.00165265],
                [+0.00083371, +0.00149466, -0.00200719, -0.00296284],
                [-0.00111959, -0.00200719, +0.00269546, +0.00397881],
                [-0.00165265, -0.00296284, +0.00397881, +0.00587320]]
    assert_array_almost_equal(expected, mmc.metric(), decimal=6)

    # Diagonal metric
    mmc = MMC(diagonal=True)
    mmc.fit(self.iris_points, [a, b, c, d])
    expected = [0, 0, 1.21045968, 1.22552608]
    assert_array_almost_equal(np.diag(expected), mmc.metric(), decimal=6)

    # Supervised Full
    mmc = MMC_Supervised()
    mmc.fit(self.iris_points, self.iris_labels)
    csep = class_separation(mmc.transform(), self.iris_labels)
    self.assertLess(csep, 0.15)

    # Supervised Diagonal
    mmc = MMC_Supervised(diagonal=True)
    mmc.fit(self.iris_points, self.iris_labels)
    csep = class_separation(mmc.transform(), self.iris_labels)
    self.assertLess(csep, 0.2)

  def __call__(self,x,v=0,t=0,calc_auxilary=True):
    distance_matrix = self.distance_matrix
    sigma,eps,m,r_tol = self.sigma,self.eps,self._masses,self.r_tol
    dx = distance_matrix(x)		# compute distance matrix THIS IS THE MAIN TIME COMPUTATION
    r = distance_matrix.radii(dx)
    bad_diagonal = isnan(r)
    if (abs(r) < r_tol).any():		# This is some exception handling for forces that get out of control
      idx = find(abs(r) < r_tol)
      err_string = "Two points within 1e-8 of eachother r = \n"
      err_string += str(r) + "\n"
      err_string += "\n r{} = \n".format(array(unravel_index(idx,r.shape)))
      for x in r.ravel()[idx]:
	err_string += "{:.12f} \n".format(float(x))
#      raise ValueError(err_string)
    K2 = (sigma/r)**6		      # Save these constants for calculating potential energy
    K1 = K2**2 
    fmatrix = 24*eps/r*(2*K1 - K2)*dx.transpose(2,0,1)   
    K1[bad_diagonal] = 0	      # Set diagonal back to zero.  
    K2[bad_diagonal] = 0	      # Set diagonal back to zero.  
    fmatrix[:,bad_diagonal] = 0	      # Set diagonal back to zero.  
    if calc_auxilary:
      self.potential_energy = sum(triu(4*eps*(K1 - K2))) # Calculate potential energy
      self.kinetic_energy = 3/2*sum(self._masses*sum(v**2,axis=1),axis=0)/2.
      N = dx.shape[0]
      d = dx.shape[2] 	     
      L = self.L[abs(self.L)>0]
      self.pressure = 0
      if N > 1:
	press_matrix = sum(fmatrix*dx.transpose(2,0,1),axis=0) # this is the "dot product"
	self.pressure = (sum(triu(press_matrix)) + N/(N-1)*self.kinetic_energy*2.)/d/prod(L)
    ans = sum(fmatrix,axis=1)/m		      # Divide by masses to get acceleration
    return ans.T