def _replace_zero_by_x_arrays(sub_arys):
    for i in range(len(sub_arys)):
        if _nx.ndim(sub_arys[i]) == 0:
            sub_arys[i] = _nx.empty(0, dtype=sub_arys[i].dtype)
        elif _nx.sometrue(_nx.equal(_nx.shape(sub_arys[i]), 0)):
            sub_arys[i] = _nx.empty(0, dtype=sub_arys[i].dtype)
    return sub_arys

def isneginf(x, y=None):
    """
    Return True where x is -infinity, and False otherwise.

    Parameters
    ----------
    x : array_like
      The input array.
    y : array_like
      A boolean array with the same shape as `x` to store the result.

    Returns
    -------
    y : ndarray
      A boolean array where y[i] = True only if x[i] = -Inf.

    See Also
    --------
    isposinf, isfinite

    Examples
    --------
    >>> np.isneginf([-np.inf, 0., np.inf])
    array([ True, False, False], dtype=bool)

    """
    if y is None:
        x = nx.asarray(x)
        y = nx.empty(x.shape, dtype=nx.bool_)
    nx.logical_and(nx.isinf(x), nx.signbit(x), y)
    return y

 def _ismethod(self, name):
     result = empty(self.shape, dtype=bool)
     res = result.flat
     for k, val in enumerate(self.flat):
         item = val.rstrip('\x00')
         res[k] = getattr(item, name)()
     return result

def mediff1d(array, to_end=None, to_begin=None):
    """Array difference with prefixed and/or appended value."""
    a = masked_array(array, copy=True)
    if a.ndim > 1:
        a.reshape((a.size,))
    (d, m, n) = (a._data, a._mask, a.size-1)
    dd = d[1:]-d[:-1]
    if m is nomask:
        dm = nomask
    else:
        dm = m[1:]-m[:-1]
    #
    if to_end is not None:
        to_end = asarray(to_end)
        nend = to_end.size
        if to_begin is not None:
            to_begin = asarray(to_begin)
            nbegin = to_begin.size
            r_data = numeric.empty((n+nend+nbegin,), dtype=a.dtype)
            r_mask = numeric.zeros((n+nend+nbegin,), dtype=bool_)
            r_data[:nbegin] = to_begin._data
            r_mask[:nbegin] = to_begin._mask
            r_data[nbegin:-nend] = dd
            r_mask[nbegin:-nend] = dm
        else:
            r_data = numeric.empty((n+nend,), dtype=a.dtype)
            r_mask = numeric.zeros((n+nend,), dtype=bool_)
            r_data[:-nend] = dd
            r_mask[:-nend] = dm
        r_data[-nend:] = to_end._data
        r_mask[-nend:] = to_end._mask
    #
    elif to_begin is not None:
        to_begin = asarray(to_begin)
        nbegin = to_begin.size
        r_data = numeric.empty((n+nbegin,), dtype=a.dtype)
        r_mask = numeric.zeros((n+nbegin,), dtype=bool_)
        r_data[:nbegin] = to_begin._data
        r_mask[:nbegin] = to_begin._mask
        r_data[nbegin:] = dd
        r_mask[nbegin:] = dm
    #
    else:
        r_data = dd
        r_mask = dm
    return masked_array(r_data, mask=r_mask)

def isneginf(x, y=None):
    """Return a boolean array y with y[i] True for x[i] = -Inf.

    If y is an array, the result replaces the contents of y.
    """
    if y is None:
        x = asarray(x)
        y = empty(x.shape, dtype=nx.bool_)
    umath.logical_and(isinf(x), signbit(x), y)
    return y

def piecewise(x, condlist, funclist, *args, **kw):
    """Return a piecewise-defined function.

    x is the domain

    condlist is a list of boolean arrays or a single boolean array
      The length of the condition list must be n2 or n2-1 where n2
      is the length of the function list.  If len(condlist)==n2-1, then
      an 'otherwise' condition is formed by |'ing all the conditions
      and inverting.

    funclist is a list of functions to call of length (n2).
      Each function should return an array output for an array input
      Each function can take (the same set) of extra arguments and
      keyword arguments which are passed in after the function list.
      A constant may be used in funclist for a function that returns a
      constant (e.g. val  and lambda x: val are equivalent in a funclist).

    The output is the same shape and type as x and is found by
      calling the functions on the appropriate portions of x.

    Note: This is similar to choose or select, except
          the the functions are only evaluated on elements of x
          that satisfy the corresponding condition.

    The result is
           |--
           |  f1(x)  for condition1
     y = --|  f2(x)  for condition2
           |   ...
           |  fn(x)  for conditionn
           |--

    """
    x = asanyarray(x)
    n2 = len(funclist)
    if not isinstance(condlist, type([])):
        condlist = [condlist]
    n = len(condlist)
    if n == n2-1:  # compute the "otherwise" condition.
        totlist = condlist[0]
        for k in range(1, n):
            totlist |= condlist[k]
        condlist.append(~totlist)
        n += 1
    if (n != n2):
        raise ValueError, "function list and condition list must be the same"
    y = empty(x.shape, x.dtype)
    for k in range(n):
        item = funclist[k]
        if not callable(item):
            y[condlist[k]] = item
        else:
            y[condlist[k]] = item(x[condlist[k]], *args, **kw)
    return y

 def _typedmethod(self, name, myiter, dtype):
     result = empty(myiter.shape, dtype=dtype)
     res = result.flat
     for k, val in enumerate(myiter):
         newval = []
         for chk in val[1:]:
             if not chk or (chk.dtype is object_ and chk.item() is None):
                 break
             newval.append(chk)
         this_str = val[0].rstrip('\x00')
         newitem = getattr(this_str,name)(*newval)
         res[k] = newitem
     return result

 def _calculatePathLength(self, save_xd):
   '''Calculates the cumulative Euclidian d-dimensional path length of the piecewise linear curve defined by a series of points save_xd
   TODO - factor this out into an 'lpcPath'-type class
   Parameters
   ----------
   save_xd : 2-dim (n*m) numpy.array of floats containing coordinates of n, ordered, m-dimensional feature points defining a
   piecewise linear curve with n-1 segments
   
   Returns
   ------- 
   lamb : 1-dim array with n ordered entries, defining the cumulative sum of segment lengths. lamb[0] = 0.
   ''' 
   it = len(save_xd)
   lamb = empty(it)
   for i in range(it):
     if i==0:
       lamb[0] = 0
     else:
       lamb[i] = lamb[i-1] + sqrt(sum((save_xd[i] - save_xd[i-1])**2))
   return lamb

def delete(arr, obj, axis=None):
    """Return a new array with sub-arrays along an axis deleted.

    Return a new array with the sub-arrays (i.e. rows or columns)
    deleted along the given axis as specified by obj

    obj may be a slice_object (s_[3:5:2]) or an integer
    or an array of integers indicated which sub-arrays to
    remove.

    If axis is None, then ravel the array first.

    Example:
    >>> arr = [[3,4,5],
    ...       [1,2,3],
    ...       [6,7,8]]

    >>> delete(arr, 1, 1)
    array([[3, 5],
           [1, 3],
           [6, 8]])
    >>> delete(arr, 1, 0)
    array([[3, 4, 5],
           [6, 7, 8]])
    """
    wrap = None
    if type(arr) is not ndarray:
        try:
            wrap = arr.__array_wrap__
        except AttributeError:
            pass


    arr = asarray(arr)
    ndim = arr.ndim
    if axis is None:
        if ndim != 1:
            arr = arr.ravel()
        ndim = arr.ndim;
        axis = ndim-1;
    if ndim == 0:
        if wrap:
            return wrap(arr)
        else:
            return arr.copy()
    slobj = [slice(None)]*ndim
    N = arr.shape[axis]
    newshape = list(arr.shape)
    if isinstance(obj, (int, long, integer)):
        if (obj < 0): obj += N
        if (obj < 0 or obj >=N):
            raise ValueError, "invalid entry"
        newshape[axis]-=1;
        new = empty(newshape, arr.dtype, arr.flags.fnc)
        slobj[axis] = slice(None, obj)
        new[slobj] = arr[slobj]
        slobj[axis] = slice(obj,None)
        slobj2 = [slice(None)]*ndim
        slobj2[axis] = slice(obj+1,None)
        new[slobj] = arr[slobj2]
    elif isinstance(obj, slice):
        start, stop, step = obj.indices(N)
        numtodel = len(xrange(start, stop, step))
        if numtodel <= 0:
            if wrap:
                return wrap(new)
            else:
                return arr.copy()
        newshape[axis] -= numtodel
        new = empty(newshape, arr.dtype, arr.flags.fnc)
        # copy initial chunk
        if start == 0:
            pass
        else:
            slobj[axis] = slice(None, start)
            new[slobj] = arr[slobj]
        # copy end chunck
        if stop == N:
            pass
        else:
            slobj[axis] = slice(stop-numtodel,None)
            slobj2 = [slice(None)]*ndim
            slobj2[axis] = slice(stop, None)
            new[slobj] = arr[slobj2]
        # copy middle pieces
        if step == 1:
            pass
        else:  # use array indexing.
            obj = arange(start, stop, step, dtype=intp)
            all = arange(start, stop, dtype=intp)
            obj = setdiff1d(all, obj)
            slobj[axis] = slice(start, stop-numtodel)
            slobj2 = [slice(None)]*ndim
            slobj2[axis] = obj
            new[slobj] = arr[slobj2]
    else: # default behavior
        obj = array(obj, dtype=intp, copy=0, ndmin=1)
        all = arange(N, dtype=intp)
        obj = setdiff1d(all, obj)
        slobj[axis] = obj
#.........这里部分代码省略.........

def insert(arr, obj, values, axis=None):
    """Return a new array with values inserted along the given axis
    before the given indices

    If axis is None, then ravel the array first.

    The obj argument can be an integer, a slice, or a sequence of
    integers.

    Example:
    >>> a = array([[1,2,3],
    ...            [4,5,6],
    ...            [7,8,9]])

    >>> insert(a, [1,2], [[4],[5]], axis=0)
    array([[1, 2, 3],
           [4, 4, 4],
           [4, 5, 6],
           [5, 5, 5],
           [7, 8, 9]])
    """
    wrap = None
    if type(arr) is not ndarray:
        try:
            wrap = arr.__array_wrap__
        except AttributeError:
            pass

    arr = asarray(arr)
    ndim = arr.ndim
    if axis is None:
        if ndim != 1:
            arr = arr.ravel()
        ndim = arr.ndim
        axis = ndim-1
    if (ndim == 0):
        arr = arr.copy()
        arr[...] = values
        if wrap:
            return wrap(arr)
        else:
            return arr
    slobj = [slice(None)]*ndim
    N = arr.shape[axis]
    newshape = list(arr.shape)
    if isinstance(obj, (int, long, integer)):
        if (obj < 0): obj += N
        if obj < 0 or obj > N:
            raise ValueError, "index (%d) out of range (0<=index<=%d) "\
                  "in dimension %d" % (obj, N, axis)
        newshape[axis] += 1;
        new = empty(newshape, arr.dtype, arr.flags.fnc)
        slobj[axis] = slice(None, obj)
        new[slobj] = arr[slobj]
        slobj[axis] = obj
        new[slobj] = values
        slobj[axis] = slice(obj+1,None)
        slobj2 = [slice(None)]*ndim
        slobj2[axis] = slice(obj,None)
        new[slobj] = arr[slobj2]
        if wrap:
            return wrap(new)
        return new

    elif isinstance(obj, slice):
        # turn it into a range object
        obj = arange(*obj.indices(N),**{'dtype':intp})

    # get two sets of indices
    #  one is the indices which will hold the new stuff
    #  two is the indices where arr will be copied over

    obj = asarray(obj, dtype=intp)
    numnew = len(obj)
    index1 = obj + arange(numnew)
    index2 = setdiff1d(arange(numnew+N),index1)
    newshape[axis] += numnew
    new = empty(newshape, arr.dtype, arr.flags.fnc)
    slobj2 = [slice(None)]*ndim
    slobj[axis] = index1
    slobj2[axis] = index2
    new[slobj] = values
    new[slobj2] = arr

    if wrap:
        return wrap(new)
    return new

def vander(x, N=None, increasing=False):
    """
    Generate a Vandermonde matrix.

    The columns of the output matrix are powers of the input vector. The
    order of the powers is determined by the `increasing` boolean argument.
    Specifically, when `increasing` is False, the `i`-th output column is
    the input vector raised element-wise to the power of ``N - i - 1``. Such
    a matrix with a geometric progression in each row is named for Alexandre-
    Theophile Vandermonde.

    Parameters
    ----------
    x : array_like
        1-D input array.
    N : int, optional
        Number of columns in the output.  If `N` is not specified, a square
        array is returned (``N = len(x)``).
    increasing : bool, optional
        Order of the powers of the columns.  If True, the powers increase
        from left to right, if False (the default) they are reversed.

        .. versionadded:: 1.9.0

    Returns
    -------
    out : ndarray
        Vandermonde matrix.  If `increasing` is False, the first column is
        ``x^(N-1)``, the second ``x^(N-2)`` and so forth. If `increasing` is
        True, the columns are ``x^0, x^1, ..., x^(N-1)``.

    See Also
    --------
    polynomial.polynomial.polyvander

    Examples
    --------
    >>> x = np.array([1, 2, 3, 5])
    >>> N = 3
    >>> np.vander(x, N)
    array([[ 1,  1,  1],
           [ 4,  2,  1],
           [ 9,  3,  1],
           [25,  5,  1]])

    >>> np.column_stack([x**(N-1-i) for i in range(N)])
    array([[ 1,  1,  1],
           [ 4,  2,  1],
           [ 9,  3,  1],
           [25,  5,  1]])

    >>> x = np.array([1, 2, 3, 5])
    >>> np.vander(x)
    array([[  1,   1,   1,   1],
           [  8,   4,   2,   1],
           [ 27,   9,   3,   1],
           [125,  25,   5,   1]])
    >>> np.vander(x, increasing=True)
    array([[  1,   1,   1,   1],
           [  1,   2,   4,   8],
           [  1,   3,   9,  27],
           [  1,   5,  25, 125]])

    The determinant of a square Vandermonde matrix is the product
    of the differences between the values of the input vector:

    >>> np.linalg.det(np.vander(x))
    48.000000000000043
    >>> (5-3)*(5-2)*(5-1)*(3-2)*(3-1)*(2-1)
    48

    """
    x = asarray(x)
    if x.ndim != 1:
        raise ValueError("x must be a one-dimensional array or sequence.")
    if N is None:
        N = len(x)

    v = empty((len(x), N), dtype=promote_types(x.dtype, int))
    tmp = v[:, ::-1] if not increasing else v

    if N > 0:
        tmp[:, 0] = 1
    if N > 1:
        tmp[:, 1:] = x[:, None]
        multiply.accumulate(tmp[:, 1:], out=tmp[:, 1:], axis=1)

    return v

def isposinf(x, y=None):
    """
    Test element-wise for positive infinity, return result as bool array.

    Parameters
    ----------
    x : array_like
        The input array.
    y : array_like, optional
        A boolean array with the same shape as `x` to store the result.

    Returns
    -------
    y : ndarray
        A boolean array with the same dimensions as the input.
        If second argument is not supplied then a boolean array is returned
        with values True where the corresponding element of the input is
        positive infinity and values False where the element of the input is
        not positive infinity.

        If a second argument is supplied the result is stored there. If the
        type of that array is a numeric type the result is represented as zeros
        and ones, if the type is boolean then as False and True.
        The return value `y` is then a reference to that array.

    See Also
    --------
    isinf, isneginf, isfinite, isnan

    Notes
    -----
    Numpy uses the IEEE Standard for Binary Floating-Point for Arithmetic
    (IEEE 754).

    Errors result if the second argument is also supplied when `x` is a
    scalar input, or if first and second arguments have different shapes.

    Examples
    --------
    >>> np.isposinf(np.PINF)
    array(True, dtype=bool)
    >>> np.isposinf(np.inf)
    array(True, dtype=bool)
    >>> np.isposinf(np.NINF)
    array(False, dtype=bool)
    >>> np.isposinf([-np.inf, 0., np.inf])
    array([False, False,  True], dtype=bool)

    >>> x = np.array([-np.inf, 0., np.inf])
    >>> y = np.array([2, 2, 2])
    >>> np.isposinf(x, y)
    array([0, 0, 1])
    >>> y
    array([0, 0, 1])

    """
    if y is None:
        x = nx.asarray(x)
        y = nx.empty(x.shape, dtype=nx.bool_)
    nx.logical_and(nx.isinf(x), ~nx.signbit(x), y)
    return y

 def test_empty_method(self):
     a = empty((2,3))
     self.assertEqual(a.ndim, 2)

def diag(v, k=0):
    """
    Extract a diagonal or construct a diagonal array.

    Parameters
    ----------
    v : array_like
        If `v` is a 2-D array, return a copy of its `k`-th diagonal.
        If `v` is a 1-D array, return a 2-D array with `v` on the `k`-th
        diagonal.
    k : int, optional
        Diagonal in question. The default is 0. Use `k>0` for diagonals
        above the main diagonal, and `k<0` for diagonals below the main
        diagonal.

    Returns
    -------
    out : ndarray
        The extracted diagonal or constructed diagonal array.

    See Also
    --------
    diagonal : Return specified diagonals.
    diagflat : Create a 2-D array with the flattened input as a diagonal.
    trace : Sum along diagonals.
    triu : Upper triangle of an array.
    tril : Lower triange of an array.

    Examples
    --------
    >>> x = np.arange(9).reshape((3,3))
    >>> x
    array([[0, 1, 2],
           [3, 4, 5],
           [6, 7, 8]])

    >>> np.diag(x)
    array([0, 4, 8])
    >>> np.diag(x, k=1)
    array([1, 5])
    >>> np.diag(x, k=-1)
    array([3, 7])

    >>> np.diag(np.diag(x))
    array([[0, 0, 0],
           [0, 4, 0],
           [0, 0, 8]])

    """
    v = asarray(v)
    s = v.shape
    if len(s) == 1:
        n = s[0]+abs(k)
        res = zeros((n,n), v.dtype)
        if k >= 0:
            i = k
        else:
            i = (-k) * n
        res[:n-k].flat[i::n+1] = v
        return res
    elif len(s) == 2:
        if k >= s[1]:
            return empty(0, dtype=v.dtype)
        if v.flags.f_contiguous:
            # faster slicing
            v, k, s = v.T, -k, s[::-1]
        if k >= 0:
            i = k
        else:
            i = (-k) * s[1]
        return v[:s[1]-k].flat[i::s[1]+1]
    else:
        raise ValueError("Input must be 1- or 2-d.")

def isposinf(x, y=None):
    """
    Shows which elements of the input are positive infinity.

    Returns a numpy array resulting from an element-wise test for positive
    infinity.

    Parameters
    ----------
    x : array_like
      The input array.
    y : array_like
      A boolean array with the same shape as `x` to store the result.

    Returns
    -------
    y : ndarray
      A numpy boolean array with the same dimensions as the input.
      If second argument is not supplied then a numpy boolean array is returned
      with values True where the corresponding element of the input is positive
      infinity and values False where the element of the input is not positive
      infinity.

      If second argument is supplied then an numpy integer array is returned
      with values 1 where the corresponding element of the input is positive
      positive infinity.

    See Also
    --------
    isinf : Shows which elements are negative or positive infinity.
    isneginf : Shows which elements are negative infinity.
    isnan : Shows which elements are Not a Number (NaN).
    isfinite: Shows which elements are not: Not a number, positive and
             negative infinity

    Notes
    -----
    Numpy uses the IEEE Standard for Binary Floating-Point for Arithmetic
    (IEEE 754). This means that Not a Number is not equivalent to infinity.
    Also that positive infinity is not equivalent to negative infinity. But
    infinity is equivalent to positive infinity.

    Errors result if second argument is also supplied with scalar input or
    if first and second arguments have different shapes.

    Numpy's definitions for positive infinity (PINF) and negative infinity
    (NINF) may be change in the future versions.


    Examples
    --------
    >>> np.isposinf(np.PINF)
    array(True, dtype=bool)
    >>> np.isposinf(np.inf)
    array(True, dtype=bool)
    >>> np.isposinf(np.NINF)
    array(False, dtype=bool)
    >>> np.isposinf([-np.inf, 0., np.inf])
    array([False, False,  True], dtype=bool)
    >>> x=np.array([-np.inf, 0., np.inf])
    >>> y=np.array([2,2,2])
    >>> np.isposinf(x,y)
    array([1, 0, 0])
    >>> y
    array([1, 0, 0])

    """
    if y is None:
        x = nx.asarray(x)
        y = nx.empty(x.shape, dtype=nx.bool_)
    nx.logical_and(nx.isinf(x), ~nx.signbit(x), y)
    return y

 def _followxSingleDirection(  self, 
                               x, 
                               direction = Direction.FORWARD,
                               forward_curve = None,
                               last_eigenvector = None, 
                               weights = 1.):
   '''Generates a partial lpc curve dictionary from the start point, x.
   Arguments
   ---------
   x : 1-dim, length m, numpy.array of floats, start point for the algorithm when m is dimension of feature space
   
   direction :  bool, proceeds in Direction.FORWARD or Direction.BACKWARD from this point (just sets sign for first eigenvalue) 
   
   forward_curve : dictionary as returned by this function, is used to detect crossing of the curve under construction with a
       previously constructed curve
       
   last_eigenvector : 1-dim, length m, numpy.array of floats, a unit vector that defines the initial direction, relative to
       which the first eigenvector is biased and initial cos_neu_neu is calculated  
       
   weights : 1-dim, length n numpy.array of observation weights (can also be used to exclude
       individual observations from the computation by setting their weight to zero.),
       where n is the number of feature points 
   '''
   x0 = copy(x)
   N = self.Xi.shape[0]
   d = self.Xi.shape[1]
   it = self._lpcParameters['it']
   h = array(self._lpcParameters['h'])
   t0 = self._lpcParameters['t0']
   rho0 = self._lpcParameters['rho0']
   
   save_xd = empty((it,d))
   eigen_vecd = empty((it,d))
   c0 = ones(it)
   cos_alt_neu = ones(it)
   cos_neu_neu = ones(it)    
   lamb = empty(it) #NOTE this is named 'lambda' in the original R code
   rho = zeros(it)
   high_rho_points = empty((0,d))    
   count_points = 0
   
   for i in range(it):
     kernel_weights = self._kernd(self.Xi, x0, c0[i]*h) * weights
     mu_x = average(self.Xi, axis = 0, weights = kernel_weights)
     sum_weights = sum(kernel_weights)
     mean_sub = self.Xi - mu_x 
     cov_x = dot( dot(transpose(mean_sub), numpy.diag(kernel_weights)), mean_sub) / sum_weights 
     #assert (abs(cov_x.transpose() - cov_x)/abs(cov_x.transpose() + cov_x) < 1e-6).all(), 'Covariance matrix not symmetric, \n cov_x = {0}, mean_sub = {1}'.format(cov_x, mean_sub)
     save_xd[i] = mu_x #save first point of the branch
     count_points += 1
     
     #calculate path length
     if i==0:
       lamb[0] = 0
     else:
       lamb[i] = lamb[i-1] + sqrt(sum((mu_x - save_xd[i-1])**2))
     
     #calculate eigenvalues/vectors
     #(sorted_eigen_cov is a list of tuples containing eigenvalue and associated eigenvector, sorted descending by eigenvalue)
     eigen_cov = eigh(cov_x)
     sorted_eigen_cov = zip(eigen_cov[0],map(ravel,vsplit(eigen_cov[1].transpose(),len(eigen_cov[1]))))
     sorted_eigen_cov.sort(key = lambda elt: elt[0], reverse = True)   
     eigen_norm = sqrt(sum(sorted_eigen_cov[0][1]**2))
     eigen_vecd[i] = direction * sorted_eigen_cov[0][1] / eigen_norm  #Unit eigenvector corresponding to largest eigenvalue
     
     #rho parameters
     rho[i] = sorted_eigen_cov[1][0] / sorted_eigen_cov[0][0] #Ratio of two largest eigenvalues
     if i != 0 and rho[i] > rho0 and rho[i-1] <= rho0:
       high_rho_points = vstack((high_rho_points, x0))
     
     #angle between successive eigenvectors
     if i==0 and last_eigenvector is not None:
       cos_alt_neu[i] = direction * dot(last_eigenvector, eigen_vecd[i])
     if i > 0:
       cos_alt_neu[i] = dot(eigen_vecd[i], eigen_vecd[i-1])
     
     #signum flipping
     if cos_alt_neu[i] < 0:
       eigen_vecd[i] = -eigen_vecd[i]
       cos_neu_neu[i] = -cos_alt_neu[i]
     else:
       cos_neu_neu[i] = cos_alt_neu[i]
    
     #angle penalization
     pen = self._lpcParameters['pen']
     if pen > 0:
       if i == 0 and last_eigenvector is not None:
         a = abs(cos_alt_neu[i])**pen
         eigen_vecd[i] = a * eigen_vecd[i] + (1-a) * last_eigenvector
       if i > 0:
         a = abs(cos_alt_neu[i])**pen
         eigen_vecd[i] = a * eigen_vecd[i] + (1-a) * eigen_vecd[i-1]
             
     #check curve termination criteria
     if i not in (0, it-1):
       #crossing
       cross = self._lpcParameters['cross']
       if forward_curve is None:
         full_curve_points = save_xd[0:i+1]
       else:
#.........这里部分代码省略.........

def histogramdd(sample, bins=10, range=None, normed=False, weights=None):
    """histogramdd(sample, bins=10, range=None, normed=False, weights=None)

    Return the N-dimensional histogram of the sample.

    Parameters:

        sample : sequence or array
            A sequence containing N arrays or an NxM array. Input data.

        bins : sequence or scalar
            A sequence of edge arrays, a sequence of bin counts, or a scalar
            which is the bin count for all dimensions. Default is 10.

        range : sequence
            A sequence of lower and upper bin edges. Default is [min, max].

        normed : boolean
            If False, return the number of samples in each bin, if True,
            returns the density.

        weights : array
            Array of weights.  The weights are normed only if normed is True.
            Should the sum of the weights not equal N, the total bin count will
            not be equal to the number of samples.

    Returns:

        hist : array
            Histogram array.

        edges : list
            List of arrays defining the lower bin edges.

    SeeAlso:

        histogram

    Example

        >>> x = random.randn(100,3)
        >>> hist3d, edges = histogramdd(x, bins = (5, 6, 7))

    """

    try:
        # Sample is an ND-array.
        N, D = sample.shape
    except (AttributeError, ValueError):
        # Sample is a sequence of 1D arrays.
        sample = atleast_2d(sample).T
        N, D = sample.shape

    nbin = empty(D, int)
    edges = D*[None]
    dedges = D*[None]
    if weights is not None:
        weights = asarray(weights)

    try:
        M = len(bins)
        if M != D:
            raise AttributeError, 'The dimension of bins must be a equal to the dimension of the sample x.'
    except TypeError:
        bins = D*[bins]

    # Select range for each dimension
    # Used only if number of bins is given.
    if range is None:
        smin = atleast_1d(array(sample.min(0), float))
        smax = atleast_1d(array(sample.max(0), float))
    else:
        smin = zeros(D)
        smax = zeros(D)
        for i in arange(D):
            smin[i], smax[i] = range[i]

    # Make sure the bins have a finite width.
    for i in arange(len(smin)):
        if smin[i] == smax[i]:
            smin[i] = smin[i] - .5
            smax[i] = smax[i] + .5

    # Create edge arrays
    for i in arange(D):
        if isscalar(bins[i]):
            nbin[i] = bins[i] + 2 # +2 for outlier bins
            edges[i] = linspace(smin[i], smax[i], nbin[i]-1)
        else:
            edges[i] = asarray(bins[i], float)
            nbin[i] = len(edges[i])+1  # +1 for outlier bins
        dedges[i] = diff(edges[i])

    nbin =  asarray(nbin)

    # Compute the bin number each sample falls into.
    Ncount = {}
    for i in arange(D):
        Ncount[i] = digitize(sample[:,i], edges[i])

#.........这里部分代码省略.........

def masked_all(shape, dtype=float_):
    """Returns an empty masked array of the given shape and dtype,
    where all the data are masked."""
    a = masked_array(numeric.empty(shape, dtype),
                     mask=numeric.ones(shape, bool_))
    return a