def intercept(self, ray):
        """Solves for intersection point of surface and a ray or Beam
    
        Args:
            ray: Ray or Beam object
                It must be in the same coordinate space as the surface object.
            
        Returns:
            s: value of s [meters] which intercepts along norm, otherwise an
            empty tuple (for no intersection).
        
        Examples:
            Accepts all point and point-derived object inputs, though all data 
            is stored as a python object.

            Generate an y direction Ray in cartesian coords using a Vec from (0,0,1)::
            
                    cen = geometry.Center(flag=True)
                    ydir = geometry.Vecx((0,1,0))
                    zpt = geometry.Point((0,0,1),cen)

        """


        # Proceedure will be to generate 
        if self._origin is ray._origin:
            try:
                rcopy = ray.copy()
                rcopy.redefine(self)
                
                intersect = _beam.interceptCyl(scipy.atleast_2d(rcopy.x()[:,-1]), 
                                               scipy.atleast_2d(rcopy.norm.unit), 
                                               scipy.array([self.sagi.s,self.sagi.s]),
                                               scipy.array([-self.norm.s,self.norm.s])) + rcopy.norm.s[-1]
                
                if not scipy.isfinite(intersect):
                    #relies on r1 using arctan2 so that it sets the branch cut properly (-pi,pi]
                    return None
                elif self.edgetest(intersect, (rcopy(intersect)).r1()):
                        return intersect
                else:
                    rcopy.norm.s[-1] = intersect
                    intersect = _beam.interceptCyl(scipy.atleast_2d(rcopy.x()[:,-1]), 
                                                   scipy.atleast_2d(rcopy.norm.unit), 
                                                   scipy.array([self.sagi.s,self.sagi.s]),
                                                   scipy.array([-self.norm.s,self.norm.s])) + rcopy.norm.s[-1]
                    if not scipy.isfinite(intersect):
                        #relies on r1 using arctan2 so that it sets the branch cut properly (-pi,pi]
                        return None
                    elif self.edgetest(intersect, (rcopy(intersect)).r1()):
                        return None
                    else:
                        return None

            except AttributeError:
                raise ValueError('not a surface object')
        else:           
            raise ValueError('not in same coordinate system, use redefine and try again')

 def __call__(self, X, n):
     n = scipy.atleast_2d(scipy.asarray(n, dtype=int))
     X = scipy.atleast_2d(scipy.asarray(X))
     n_unique = unique_rows(n)
     mu = scipy.zeros(X.shape[0])
     for nn in n_unique:
         idxs = (n == nn).all(axis=1)
         mu[idxs] = self.fun(X[idxs, :], nn, *self.params)
     
     return mu

	def topterms(self,n_terms=10):
		""" This function is given. """
		vec = sp.atleast_2d(sp.arange(0,self.n_words))
		topics = []
		for k in xrange(self.n_topics):
			probs = sp.atleast_2d(self._phi[k,:])
			mat = sp.append(probs,vec,0)
			sind = sp.array([mat[:,i] for i in sp.argsort(mat[0])]).T
			topics.append([self.vocab[int(sind[1,self.n_words - 1 - i])] for i in xrange(n_terms)])
		return topics

def worker_quality(predictions, num_classes):
    predictions = sp.atleast_2d(predictions)
    num_workers, num_objects = predictions.shape

    error_rates = sp.zeros((num_workers, num_classes, num_classes))
    diy, diz = sp.diag_indices(num_classes)
    error_rates[:, diy, diz] = 1

    while True:
        # E step
        new_predictions = sp.zeros((num_objects, num_classes))
        for i in xrange(num_objects):
            individual_predictions = predictions[:, i]
            individual_error_rates = error_rates[range(num_workers), individual_predictions, individual_predictions]
            new_predictions[i, :] = sp.bincount(individual_predictions, individual_error_rates, minlength=num_classes)

        correct_labels = sp.argmax(new_predictions, axis=1)
        count_per_label = sp.bincount(correct_labels)

        # M step
        new_error_rates = sp.zeros((num_workers, num_classes, num_classes))
        for i, label in enumerate(correct_labels):
            new_error_rates[range(num_workers), label, predictions[:, i]] += 1

        for i in xrange(num_classes):
            new_error_rates[:, :, i] /= count_per_label

        diff_error_rates = sp.absolute(new_error_rates - error_rates)
        error_rates = new_error_rates

        if sp.amax(diff_error_rates) < 0.001:
            break


    # calculate the cost of each worker
    class_priors = sp.bincount(correct_labels, minlength=num_classes) / float(num_objects)
    costs = []
    for k in xrange(num_workers):
        worker_class_priors = sp.dot(sp.atleast_2d(class_priors), error_rates[k])[0] + 0.0000001

        cost = 0
        for j in xrange(num_classes):
            soft_label = error_rates[k, :, j] * class_priors / worker_class_priors[j]

            soft_label_cost = 0.0
            for i in xrange(num_classes):
                soft_label_cost += sp.sum(soft_label[i] * soft_label)
            soft_label_cost -= sp.sum(soft_label ** 2) # subtract the diagonal entries (those costs = 0)
            cost += soft_label_cost * worker_class_priors[j]

        costs.append(cost)

    return error_rates, correct_labels, costs

def vec2ten(data, nchan=4):
    """converts from templates/spikes that are concatenated across the
    channels to tensors that have an extra dim for the channels

    :type data: ndarray
    :param data: input array [templates][vars * channels]
    :type nchan: int
    :param nchan: count of channels
        Default=4
    :returns: ndarray - data converted to tensor [templates][vars][channels]
    """

    if data.ndim == 1:
        data = sp.atleast_2d(data)
    n, dim = data.shape

    if dim % nchan != 0:
        raise ValueError(
            'dim %s nchan != 0 !! dim=%s, nchan=%s' % (dim, nchan))
    tf = dim / nchan

    rval = sp.zeros((n, tf, nchan), data.dtype)

    for i in xrange(n):
        for c in xrange(nchan):
            rval[i, :, c] = data[i, c * tf:(c + 1) * tf]
    return rval

def nullspace(A, atol=1e-13, rtol=0):
    '''Compute an approximate basis for the nullspace of A.
    The algorithm used by this function is based on the singular value
    decomposition of `A`. This implementation was copied
    from the scipy cookbook: http://www.scipy.org/Cookbook/RankNullspace

    @param A: ndarray
        A should be at most 2-D.  A 1-D array with length k will be treated
        as a 2-D with shape (1, k)
    @param atol : float
        The absolute tolerance for a zero singular value.  Singular values
        smaller than `atol` are considered to be zero.
    @param rtol : float
        The relative tolerance.  Singular values less than rtol*smax are
        considered to be zero, where smax is the largest singular value.

    @note: If both `atol` and `rtol` are positive, the combined tolerance is the
    maximum of the two; that is::
        tol = max(atol, rtol * smax)
    Singular values smaller than `tol` are considered to be zero.

    @return: ns ndarray
        If `A` is an array with shape (m, k), then `ns` will be an array
        with shape (k, n), where n is the estimated dimension of the
        nullspace of `A`.  The columns of `ns` are a basis for the
        nullspace; each element in numpy.dot(A, ns) will be approximately
        zero.
    '''

    A = sp.atleast_2d(A)
    _u, s, vh = LA.svd(A)
    tol = max(atol, rtol * s[0])
    nnz = (s >= tol).sum()
    ns = vh[nnz:].conj().T
    return ns

def sinc_interp1d(x, s, r):
    """Interpolates `x`, sampled at times `s`
    Output `y` is sampled at times `r`

    inspired from from Matlab:
    http://phaseportrait.blogspot.com/2008/06/sinc-interpolation-in-matlab.html

    :param ndarray x: input data time series
    :param ndarray s: input sampling time series (regular sample interval)
    :param ndarray r: output sampling time series
    :return ndarray: output data time series (regular sample interval)
    """

    # init
    s = sp.asarray(s)
    r = sp.asarray(r)
    x = sp.asarray(x)
    if x.ndim == 1:
        x = sp.atleast_2d(x)
    else:
        if x.shape[0] == len(s):
            x = x.T
        else:
            if x.shape[1] != s.shape[0]:
                raise ValueError('x and s must be same temporal extend')
    if sp.allclose(s, r):
        return x.T
    T = s[1] - s[0]

    # resample
    sincM = sp.tile(r, (len(s), 1)) - sp.tile(s[:, sp.newaxis], (1, len(r)))
    return sp.vstack([sp.dot(xx, sp.sinc(sincM / T)) for xx in x]).T

    def reconstruct(self, X):
        n_features = sp.atleast_2d(X).shape[1]
        latent = sp.dot(self.inv_M, sp.dot(self.weight.T, (X - self.predict_mean).T))
        eps = sprd.multivariate_normal(sp.zeros(n_features), self.sigma2 * sp.eye(n_features))
        recons = sp.dot(self.weight, latent) + self.predict_mean + eps

        return recons

    def phigprov(self, Pp, Pg, theta):
        """ Calculate transition probabilities

        Parameters
        ------------
        Pp : ndarray, shape (n, k)
             Conditional choice probabilities for provinces
        Pg : ndarray, shape (n, 2 k)
             Conditional choice probabilities for the government
        theta : ndarray, shape (5, )
             Parameters

        Returns
        ---------
        V : ndarray
            Observable state values

        Notes
        -----------

        Takes conditional choice probabilities :math:`P` and :math:`\theta`
        as an input and returns values :math:`V^P`.
        This is the mapping :math:`\Phi` in part (b) of the assignment.

        This is a wrapper for the matlab function **Phigprov**.
        
        """
        theta = sp.atleast_2d(theta)
        return pytave.feval(1, "Phigprov", Pp, Pg, theta, self.model())[0]

    def plot_filter_set(self, ph=None, show=False):
        """plot the filter set in a waveform plot"""

        # get plotting tools
        try:
            from spikeplot import waveforms
        except ImportError:
            return None

        # checks
        if self.nf == 0:
            warnings.warn("skipping plot, no active units!")
            return None

        # init
        units = {}
        for k in self._idx_active_set:
            units[k] = sp.atleast_2d(self.bank[k].f_conc)

        return waveforms(
            units,
            tf=self._tf,
            plot_separate=True,
            plot_mean=False,
            plot_single_waveforms=False,
            plot_handle=ph,
            show=show,
        )

    def new_p(self, Pp, Pg, theta):
        """ Calculate transition probabilities

        Parameters
        --------------
        
        Pp : ndarray, shape (n, k)
             Conditional choice probabilities for provinces
        Pg : ndarray, shape (n, 2 k)
             Conditional choice probabilities for the government
        theta : ndarray, shape (5, )
             Parameters

        Returns
        ---------
        Pp : ndarray, shape (n, k)
             New conditional choice probabilities for provinces
        Pg : ndarray, shape (n, 2 k)
             New conditional choice probabilities for the government

        Notes
        -----------

        Takes conditional choice probabilities :math:`P` and :math:`\theta`
        as an input and returns new conditional choice values.
        This is the mapping :math:`\Psi` in part (c) of the assignment.

        This is a wrapper for the matlab function **NewP**.
        
        """
        theta = sp.atleast_2d(theta)
        return pytave.feval(2, "NewP", Pp, Pg, theta, self.model())

def chunk_data(data, epochs=None, invert=False):
    """returns a generator of chunks from data given epochs

    :type data: ndarray
    :param data: signal data [[samples, channels]]
    :type epochs: ndarray
    :param epochs: epoch set, positive mask
    :type invert: bool
    :param invert: invert epochs, negative mask instead of positive mask
    :returns: generator - data chunks as per :epochs:
    """

    # checks
    data = sp.asarray(data)
    if data.ndim != 2:
        data = sp.atleast_2d(data).T
    if epochs is not None:
        if epochs.ndim != 2:
            raise ValueError("epochs has to be ndim=2 like [[start,end]]")
    if invert is True and epochs is not None:
        epochs = invert_epochs(epochs, end=data.shape[0])
    if epochs is None or len(epochs) == 0:
        epochs = [[0, data.shape[0]]]

    # yield data chunks
    for ep in epochs:
        yield data[ep[0] : ep[1], :], list(ep)

    def summed_dist_matrix(self, vectors, presorted=False):
        """ Calculates the sum of all element pair distances for each
        pair of vectors.

        If :math:`(a_1, \\dots, a_n)` and :math:`(b_1, \\dots, b_m)` are the
        :math:`u`-th and :math:`v`-th vector from `vectors` and :math:`K` the
        kernel, the resulting entry in the 2D array will be :math:`D_{uv}
        = \\sum_{i=1}^{n} \\sum_{j=1}^{m} K(a_i - b_j)`.

        :param sequence vectors: A sequence of Quantity 1D to calculate the
            summed distances for each pair. The required units depend on the
            kernel. Usually it will be the inverse unit of the kernel size.
        :param bool presorted: Some optimized specializations of this function
            may need sorted vectors. Set `presorted` to `True` if you know that
            the passed vectors are already sorted to skip the sorting and thus
            increase performance.
        :rtype: Quantity 2D
        """

        D = sp.empty((len(vectors), len(vectors)))
        if len(vectors) > 0:
            might_have_units = self(vectors[0])
            if hasattr(might_have_units, 'units'):
                D = D * might_have_units.units
            else:
                D = D * pq.dimensionless

        for i, j in sp.ndindex(len(vectors), len(vectors)):
            D[i, j] = sp.sum(self(
                (vectors[i] - sp.atleast_2d(vectors[j]).T).flatten()))
        return D

def signal(signal, events=None, epochs=None, spike_trains=None,
                spike_waveforms=None):
    """ Create a plot from an AnalogSignal.

    :param AnalogSignal signal: The signal to plot.
    :param sequence events: A list of Event objects to be included in the
        plot.
    :param sequence epochs: A list of Epoch objects to be included in the
        plot.
    :param dict spike_trains: A dictionary of SpikeTrain objects to be
        included in the plot. Spikes are plotted as vertical lines.
        Indices of the dictionary (typically Unit objects) are used
        for color and legend entries.
    :param sequence spike_waveforms: A dictionary of lists of Spike objects
        to be included in the plot. Waveforms of spikes are overlaid on
        the signal. Indices of the dictionary (typically Unit objects) are
        used for color and legend entries.
    """
    # Plot title
    win_title = 'Analog Signal'
    if signal.recordingchannel:
        win_title += ' | Recording Channel: %s' % \
                     signal.recordingchannel.name
    if signal.segment:
        win_title += ' | Segment: %s' % signal.segment.name
    win = PlotDialog(toolbar=True, wintitle=win_title)

    signalarray = neo.AnalogSignalArray(sp.atleast_2d(sp.asarray(signal)).T,
        units=signal.units, sampling_rate=signal.sampling_rate)

    _plot_signal_array_on_window(win, signalarray, events, epochs,
        spike_trains, spike_waveforms, False)

 def _stop_training(self, *args, **kwargs):
     # produce data in one piece
     self.data = sp.vstack(self.data)
     # calculate energy
     self.energy = self._energy_func(self.data)
     if self.energy.ndim == 1:
         self.energy = sp.atleast_2d(self.energy).T
     self.size, self.nchan = self.energy.shape

def genCylGrid(x0,x1,x2,edges=False):

    if edges:
        for i in (x0,x1,x2):
            i = scipy.insert(i,0,2*i[1]-i[2])
            i = scipy.append(i,2*i[-1]-i[-2])
            i = (i[1:]+i[:-1])/2

    pnts = scipy.empty((x0.size, x1.size, x2.size,3))
    xin = scipy.dot(scipy.atleast_2d(x0).T, scipy.atleast_2d(scipy.cos(x1)))
    yin = scipy.dot(scipy.atleast_2d(x0).T, scipy.atleast_2d(scipy.sin(x1)))
    zee = scipy.ones(yin.shape)
    for i in range(x2.size):
        pnts[:,:,i,0] = xin
        pnts[:,:,i,1] = yin
        pnts[:,:,i,2] = x2[i]*zee
    return pnts

def _multi_norm(x, mean):
    """ Evaluate pdf of multivariate normal distribution with a mean
        at rows of x with high precision.
    """
    d = x.shape[1]
    fac = (2 * sp.pi) ** (-d / 2.0)
    y = cdist(x, sp.atleast_2d(mean), "sqeuclidean") * -0.5
    return fac * sp.exp(sp.longdouble(y))

def get_all_forces(atoms, mlmodel, grid_spacing, extfield=None, mixing=[1.,0.,0.], lammpsdata=None, do_update=False):
    # get actual forces and potential energy of configuration
    pot_energy, forces = calc_lammps(atoms, preloaded_data=lammpsdata)
    forces = [forces]

    ### ML IS HERE ###
    if not (mlmodel is None or (mixing[1] == 0 and not do_update)):
        # Accumulate the new observation in the dataset
        coarse_colvars = round_vector(atoms.colvars(), precision=grid_spacing)
        distance_from_data = sp_dist.cdist(
            sp.atleast_2d(coarse_colvars), mlmodel.X_fit_).ravel()
        # check if configuration has already occurred
        if distance_from_data.min() == 0.0:
            index = list(distance_from_data).index(0.0)
            # Learn E_min: uncomment this
            if pot_energy < mlmodel.y[index]: 
                mlmodel.y[index] = pot_energy
            else:
                do_update = False
    #             # Learn free energy: uncomment this
    #             beta = 1 / self.temp
    #             mlmodel.y[index] += - 1 / beta * sp.log(
    #                 1 + sp.exp(- beta * (pot_energy - mlmodel.y[index])))
    #             pot_energy = mlmodel.y[index]
        else:
            mlmodel.accumulate_data(coarse_colvars, pot_energy)
        if do_update:
            # update ML potential with all the data contained in it.
            mlmodel.update_fit()
        # Get ML constraint forces if the model is fitted
        if hasattr(mlmodel, 'dual_coef_') and pot_energy < 0:
            ml_forces = get_constraint_forces(atoms, mlmodel)
            # ml_forces /= sp.mean(map(LA.norm, ml_forces))
            # ml_forces *= 1 * sp.mean(map(LA.norm, forces[0]))
            forces.append(ml_forces)
        else:
            forces.append(sp.zeros(forces[0].shape))
        # X_near = Xplotgrid([atoms.phi() - 0.2, atoms.psi() - 0.2], [atoms.phi() - 0.2, atoms.psi() - 0.2], 2, 10)
        # y_near_mean = mlmodel.predict(X_near).mean()
        # if pot_energy < y_near_mean:
        #             mix = 1.2
        #         else:
        #             mix = 0.8

    # EXTERNAL FIELD IS HERE
    if not (extfield is None or mixing[2] == 0):
        colvars = round_vector(atoms.colvars(), precision=grid_spacing)
        extfield.update_cost(colvars, pot_energy) 
        extfield_forces = get_extfield_forces(atoms, extfield)
        extfield_forces /= sp.mean(map(LA.norm, extfield_forces))
        extfield_forces *= sp.mean(map(LA.norm, forces[0]))
        forces.append(extfield_forces)
    # Compose the actual and the ML forces together by mixing them accordingly
    # a [1,-1,0] mixing would result, in the perfect fitting limit, to a zero
    # mean field motion.
    forces = [m_i * f_i for m_i, f_i in zip(mixing[:len(forces)], forces)]
    f = sp.sum(forces, axis=0)
    return f, pot_energy, forces

def regular_cube_innerproduct(rcc,k):      
    """
    For a given regular_cube_complex, compute a matrix
    representing the k-form innerproduct.

    These elements are similar to Whitney forms,
    except using standard linear (bilinear,trilinear,..)
    elements for 0-forms.
    """

    N = rcc.complex_dimension()

    #standard cube is [0,0,..,0] [0,1,...,N]   
    standard_cube  = atleast_2d(array([0]*N + range(N),dtype='i'))
    standard_k_faces = standard_cube
    for i in range(N,k,-1):        
        standard_k_faces = cube_array_boundary(standard_k_faces,i)[0]

        
    k_faces_per_cube = standard_k_faces.shape[0]


    K = zeros((k_faces_per_cube,k_faces_per_cube)) #local stiffness matrix
    h = 1
    V = h**N #cube volume
    scale = V * (1/h)**2 * (1/3.0)**(N-k)
    for i,row_i in enumerate(standard_k_faces):
        for j,row_j in enumerate(standard_k_faces):
            if all(row_i[N:] == row_j[N:]):
                differences = (row_i[:N] != row_j[:N])
                differences[row_i[N:]] = 0                
                K[i,j] = scale * (1.0/2.0)**sum(differences)
            else:
                K[i,j] = 0
        

    CA = rcc[-1].cube_array[:,:N]
    num_cubes = CA.shape[0]

    k_faces  = tile(hstack((CA,zeros((CA.shape[0],k),dtype=CA.dtype))),(1,k_faces_per_cube)).reshape((-1,N+k))
    k_faces += tile(standard_k_faces,(num_cubes,1))
    
    k_face_array = rcc[k].cube_array

    face_indices = cube_array_search(k_face_array,k_faces)

    rows = face_indices.repeat(k_faces_per_cube)
    cols = face_indices.reshape((-1,k_faces_per_cube)).repeat(k_faces_per_cube,axis=0).reshape((-1,))
    data = K.reshape((1,-1)).repeat(num_cubes,axis=0).reshape((-1,))
    
    # temporary memory cost solution - eliminate zeros from COO representation
    nz_mask = data != 0.0
    rows = rows[nz_mask]
    cols = cols[nz_mask]
    data = data[nz_mask]

    shape = (len(k_face_array),len(k_face_array))
    return coo_matrix( (data,(rows,cols)), shape).tocsr()

def detect_skew(img, min_angle=-20, max_angle=20, quality='low'):
    img = sp.atleast_2d(img)
    rows, cols = img.shape
    min_min_angle = min_angle
    max_max_angle = max_angle

    if quality == 'low':
        resolution = sp.arctan2(2.0, cols) * 180.0 / sp.pi
        min_target_size = 100
        resize_order = 1
    elif quality == 'high':
        resolution = sp.arctan2(1.0, cols) * 180.0 / sp.pi
        min_target_size = 300
        resize_order = 3
    else:
        resolution = sp.arctan2(1.0, cols) * 180.0 / sp.pi
        min_target_size = 200
        resize_order = 2

    # resize the image so it's faster to work with
    min_size = min(rows, cols)
    target_size = min_target_size if min_size > min_target_size else min_size
    resize_ratio = float(target_size) / min_size
    img = imresize(img, resize_ratio)
    rows, cols = img.shape

    # pad the image and invert the colors
    img *= -1
    img += 255
    padded_img = sp.zeros((rows*2, cols*2))
    padded_img[rows//2:rows//2+rows, cols//2:cols//2+cols] = img
    img = padded_img

    # keep dividing the interval in half to achieve O(log(n))
    while True:
        current_resolution = (max_angle - min_angle) / 30.0
        best_angle = None
        best_variance = 0.0

        # rotate the image, sum the pixel values in each row for each rotation
        # then find the variance of all the sums, pick the highest variance
        for i in xrange(31):
            angle = min_angle + i * current_resolution
            rotated_img = rotate(img, angle, reshape=False, order=resize_order)
            num_black_pixels = sp.sum(rotated_img, axis=1)
            variance = sp.var(num_black_pixels)
            if variance > best_variance:
                best_angle = angle
                best_variance = variance

        if current_resolution < resolution:
            break

        # update the angle range
        min_angle = max(best_angle - current_resolution, min_min_angle)
        max_angle = min(best_angle + current_resolution, max_max_angle)

    return best_angle