def _asarray1d(arr, copy=False):
    """Ensure 1D array for one array.
    """
    if copy:
        return asarray(arr).flatten()
    else:
        return asarray(arr).ravel()

def meshgrid(x,y):
    """
    For vectors x, y with lengths Nx=len(x) and Ny=len(y), return X, Y
    where X and Y are (Ny, Nx) shaped arrays with the elements of x
    and y repeated to fill the matrix

    EG,

      [X, Y] = meshgrid([1,2,3], [4,5,6,7])

       X =
         1   2   3
         1   2   3
         1   2   3
         1   2   3


       Y =
         4   4   4
         5   5   5
         6   6   6
         7   7   7
  """
    x = asarray(x)
    y = asarray(y)
    numRows, numCols = len(y), len(x)  # yes, reversed
    x = x.reshape(1,numCols)
    X = x.repeat(numRows, axis=0)

    y = y.reshape(numRows,1)
    Y = y.repeat(numCols, axis=1)
    return X, Y

def apply_along_axis(func1d,axis,arr,*args):
    """ Execute func1d(arr[i],*args) where func1d takes 1-D arrays
        and arr is an N-d array.  i varies so as to apply the function
        along the given axis for each 1-d subarray in arr.
    """
    arr = asarray(arr)
    nd = arr.ndim
    if axis < 0:
        axis += nd
    if (axis >= nd):
        raise ValueError("axis must be less than arr.ndim; axis=%d, rank=%d."
            % (axis,nd))
    ind = [0]*(nd-1)
    i = zeros(nd,'O')
    indlist = range(nd)
    indlist.remove(axis)
    i[axis] = slice(None,None)
    outshape = asarray(arr.shape).take(indlist)
    i.put(indlist, ind)
    res = func1d(arr[tuple(i.tolist())],*args)
    #  if res is a number, then we have a smaller output array
    if isscalar(res):
        outarr = zeros(outshape,asarray(res).dtype)
        outarr[tuple(ind)] = res
        Ntot = product(outshape)
        k = 1
        while k < Ntot:
            # increment the index
            ind[-1] += 1
            n = -1
            while (ind[n] >= outshape[n]) and (n > (1-nd)):
                ind[n-1] += 1
                ind[n] = 0
                n -= 1
            i.put(indlist,ind)
            res = func1d(arr[tuple(i.tolist())],*args)
            outarr[tuple(ind)] = res
            k += 1
        return outarr
    else:
        Ntot = product(outshape)
        holdshape = outshape
        outshape = list(arr.shape)
        outshape[axis] = len(res)
        outarr = zeros(outshape,asarray(res).dtype)
        outarr[tuple(i.tolist())] = res
        k = 1
        while k < Ntot:
            # increment the index
            ind[-1] += 1
            n = -1
            while (ind[n] >= holdshape[n]) and (n > (1-nd)):
                ind[n-1] += 1
                ind[n] = 0
                n -= 1
            i.put(indlist, ind)
            res = func1d(arr[tuple(i.tolist())],*args)
            outarr[tuple(i.tolist())] = res
            k += 1
        return outarr

 def __gmmEm__(self):
     self.mean = kmeans2(self.data, self.K)[0]
     self.c = asarray([1.0/self.K]*self.K)
     self.covm = asarray([identity(self.K)]*self.K)
     self.p = ndarray((self.N,self.K),dtype='float32')
     while self.it > 0:
         self.it -=1
         self.__calculateP__()
         #self.__Estep__()
         self.__Mstep__()

def polyval(p, x):
    """
    Evaluate a polynomial at specific values.

    If p is of length N, this function returns the value:

        p[0]*(x**N-1) + p[1]*(x**N-2) + ... + p[N-2]*x + p[N-1]

    If x is a sequence then p(x) will be returned for all elements of x.
    If x is another polynomial then the composite polynomial p(x) will
    be returned.

    Parameters
    ----------
    p : {array_like, poly1d}
       1D array of polynomial coefficients from highest degree to zero or an
       instance of poly1d.
    x : {array_like, poly1d}
       A number, a 1D array of numbers, or an instance of poly1d.

    Returns
    -------
    values : {ndarray, poly1d}
       If either p or x is an instance of poly1d, then an instance of poly1d
       is returned, otherwise a 1D array is returned. In the case where x is
       a poly1d, the result is the composition of the two polynomials, i.e.,
       substitution is used.

    See Also
    --------
    poly1d: A polynomial class.

    Notes
    -----
    Horner's method is used to evaluate the polynomial. Even so, for
    polynomials of high degree the values may be inaccurate due to
    rounding errors. Use carefully.


    Examples
    --------
    >>> np.polyval([3,0,1], 5)  # 3 * 5**2 + 0 * 5**1 + 1
    76

    """
    p = NX.asarray(p)
    if isinstance(x, poly1d):
        y = 0
    else:
        x = NX.asarray(x)
        y = NX.zeros_like(x)
    for i in range(len(p)):
        y = x * y + p[i]
    return y

def _fix_real_abs_gt_1(x):
    """Convert `x` to complex if it has real components x_i with abs(x_i)>1.

    Otherwise, output is just the array version of the input (via asarray).

    Parameters
    ----------
    x : array_like

    Returns
    -------
    array

    Examples
    --------
    >>> np.lib.scimath._fix_real_abs_gt_1([0,1])
    array([0, 1])

    >>> np.lib.scimath._fix_real_abs_gt_1([0,2])
    array([ 0.+0.j,  2.+0.j])
    """
    x = asarray(x)
    if any(isreal(x) & (abs(x)>1)):
        x = _tocomplex(x)
    return x

def mintypecode(typechars,typeset='GDFgdf',default='d'):
    """ Return a minimum data type character from typeset that
    handles all typechars given

    The returned type character must be the smallest size such that
    an array of the returned type can handle the data from an array of
    type t for each t in typechars (or if typechars is an array,
    then its dtype.char).

    If the typechars does not intersect with the typeset, then default
    is returned.

    If t in typechars is not a string then t=asarray(t).dtype.char is
    applied.
    """
    typecodes = [(type(t) is type('') and t) or asarray(t).dtype.char\
                 for t in typechars]
    intersection = [t for t in typecodes if t in typeset]
    if not intersection:
        return default
    if 'F' in intersection and 'd' in intersection:
        return 'D'
    l = []
    for t in intersection:
        i = _typecodes_by_elsize.index(t)
        l.append((i,t))
    l.sort()
    return l[0][1]

def vander(x, N=None):
    """
    Generate a Van der Monde matrix.

    The columns of the output matrix are decreasing powers of the input
    vector.  Specifically, the i-th output column is the input vector to
    the power of ``N - i - 1``. Such a matrix with a geometric progression
    in each row is named Van Der Monde, or Vandermonde matrix, from
    Alexandre-Theophile Vandermonde.

    Parameters
    ----------
    x : array_like
        1-D input array.
    N : int, optional
        Order of (number of columns in) the output. If `N` is not specified,
        a square array is returned (``N = len(x)``).

    Returns
    -------
    out : ndarray
        Van der Monde matrix of order `N`.  The first column is ``x^(N-1)``,
        the second ``x^(N-2)`` and so forth.

    References
    ----------
    .. [1] Wikipedia, "Vandermonde matrix",
           http://en.wikipedia.org/wiki/Vandermonde_matrix

    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]])

    """
    x = asarray(x)
    if N is None:
        N = len(x)
    X = ones((len(x), N), x.dtype)
    for i in range(N - 1):
        X[:, i] = x ** (N - i - 1)
    return X

def diagflat(v,k=0):
    """Return a 2D array whose k'th diagonal is a flattened v and all other
    elements are zero.

    Examples
    --------
      >>> diagflat([[1,2],[3,4]]])
      array([[1, 0, 0, 0],
             [0, 2, 0, 0],
             [0, 0, 3, 0],
             [0, 0, 0, 4]])

      >>> diagflat([1,2], 1)
      array([[0, 1, 0],
             [0, 0, 2],
             [0, 0, 0]])
    """
    try:
        wrap = v.__array_wrap__
    except AttributeError:
        wrap = None
    v = asarray(v).ravel()
    s = len(v)
    n = s + abs(k)
    res = zeros((n,n), v.dtype)
    if (k>=0):
        i = arange(0,n-k)
        fi = i+k+i*n
    else:
        i = arange(0,n+k)
        fi = i+(i-k)*n
    res.flat[fi] = v
    if not wrap:
        return res
    return wrap(res)

def diag(v, k=0):
    """ returns a copy of the the k-th diagonal if v is a 2-d array
        or returns a 2-d array with v as the k-th diagonal if v is a
        1-d array.
    """
    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 = arange(0,n-k)
            fi = i+k+i*n
        else:
            i = arange(0,n+k)
            fi = i+(i-k)*n
        res.flat[fi] = v
        return res
    elif len(s)==2:
        N1,N2 = s
        if k >= 0:
            M = min(N1,N2-k)
            i = arange(0,M)
            fi = i+k+i*N2
        else:
            M = min(N1+k,N2)
            i = arange(0,M)
            fi = i + (i-k)*N2
        return v.flat[fi]
    else:
        raise ValueError, "Input must be 1- or 2-d."

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 _fix_real_lt_zero(x):
    """Convert `x` to complex if it has real, negative components.

    Otherwise, output is just the array version of the input (via asarray).

    Parameters
    ----------
    x : array_like

    Returns
    -------
    array

    Examples
    --------
    >>> np.lib.scimath._fix_real_lt_zero([1,2])
    array([1, 2])

    >>> np.lib.scimath._fix_real_lt_zero([-1,2])
    array([-1.+0.j,  2.+0.j])
    """
    x = asarray(x)
    if any(isreal(x) & (x<0)):
        x = _tocomplex(x)
    return x

def polyint(p, m=1, k=None):
    """Return the mth analytical integral of the polynomial p.

    If k is None, then zero-valued constants of integration are used.
    otherwise, k should be a list of length m (or a scalar if m=1) to
    represent the constants of integration to use for each integration
    (starting with k[0])
    """
    m = int(m)
    if m < 0:
        raise ValueError, "Order of integral must be positive (see polyder)"
    if k is None:
        k = NX.zeros(m, float)
    k = atleast_1d(k)
    if len(k) == 1 and m > 1:
        k = k[0]*NX.ones(m, float)
    if len(k) < m:
        raise ValueError, \
              "k must be a scalar or a rank-1 array of length 1 or >m."
    if m == 0:
        return p
    else:
        truepoly = isinstance(p, poly1d)
        p = NX.asarray(p)
        y = NX.zeros(len(p)+1, float)
        y[:-1] = p*1.0/NX.arange(len(p), 0, -1)
        y[-1] = k[0]
        val = polyint(y, m-1, k=k[1:])
        if truepoly:
            val = poly1d(val)
        return val

def _fix_int_lt_zero(x):
    """Convert `x` to double if it has real, negative components.

    Otherwise, output is just the array version of the input (via asarray).

    Parameters
    ----------
    x : array_like

    Returns
    -------
    array

    Examples
    --------
    >>> _fix_int_lt_zero([1,2])
    array([1, 2])

    >>> _fix_int_lt_zero([-1,2])
    array([-1.,  2.])
    """
    x = asarray(x)
    if any(isreal(x) & (x < 0)):
        x = x * 1.0
    return x

def isrealobj(x):
    """
    Return True if x is a not complex type or an array of complex numbers.

    The type of the input is checked, not the value. So even if the input
    has an imaginary part equal to zero, `isrealobj` evaluates to False
    if the data type is complex.

    Parameters
    ----------
    x : any
        The input can be of any type and shape.

    Returns
    -------
    y : bool
        The return value, False if `x` is of a complex type.

    See Also
    --------
    iscomplexobj, isreal

    Examples
    --------
    >>> np.isrealobj(1)
    True
    >>> np.isrealobj(1+0j)
    False
    >>> np.isrealobj([3, 1+0j, True])
    False

    """
    return not issubclass(asarray(x).dtype.type, _nx.complexfloating)

def asfarray(a, dtype=_nx.float_):
    """
    Return an array converted to a float type.

    Parameters
    ----------
    a : array_like
        The input array.
    dtype : str or dtype object, optional
        Float type code to coerce input array `a`.  If `dtype` is one of the
        'int' dtypes, it is replaced with float64.

    Returns
    -------
    out : ndarray
        The input `a` as a float ndarray.

    Examples
    --------
    >>> np.asfarray([2, 3])
    array([2.,  3.])
    >>> np.asfarray([2, 3], dtype='float')
    array([2.,  3.])
    >>> np.asfarray([2, 3], dtype='int8')
    array([2.,  3.])

    """
    if not _nx.issubdtype(dtype, _nx.inexact):
        dtype = _nx.float_
    return asarray(a, dtype=dtype)

def ix_(*args):
    """ Construct an open mesh from multiple sequences.

    This function takes n 1-d sequences and returns n outputs with n
    dimensions each such that the shape is 1 in all but one dimension and
    the dimension with the non-unit shape value cycles through all n
    dimensions.

    Using ix_() one can quickly construct index arrays that will index
    the cross product.

    a[ix_([1,3,7],[2,5,8])]  returns the array

    a[1,2]  a[1,5]  a[1,8]
    a[3,2]  a[3,5]  a[3,8]
    a[7,2]  a[7,5]  a[7,8]
    """
    out = []
    nd = len(args)
    baseshape = [1]*nd
    for k in range(nd):
        new = _nx.asarray(args[k])
        if (new.ndim != 1):
            raise ValueError, "Cross index must be 1 dimensional"
        if issubclass(new.dtype.type, _nx.bool_):
            new = new.nonzero()[0]
        baseshape[k] = len(new)
        new = new.reshape(tuple(baseshape))
        out.append(new)
        baseshape[k] = 1
    return tuple(out)

 def __getitem__(self, indx):
     reset_full = True
     if isinstance(indx, Date):
         indx = self.find_dates(indx)
         reset_full = False
     elif numeric.asarray(indx).dtype.kind == 'O':
         try:
             indx = self.find_dates(indx)
         except AttributeError:
             pass
     r = ndarray.__getitem__(self, indx)
     if isinstance(r, (generic, int)):
         return Date(self.freq, value=r)
     elif hasattr(r, 'size') and r.size == 1:
         # need to check if it has a size attribute for situations
         # like when the datearray is the data for a maskedarray
         # or some other subclass of ndarray with wierd getitem
         # behaviour
         return Date(self.freq, value=r.item())
     else:
         if hasattr(r, '_cachedinfo'):
             _cache = r._cachedinfo
             _cache.update(dict([(k,_cache[k][indx])
                                 for k in ('toobj', 'tostr', 'toord')
                                 if _cache[k] is not None]))
             _cache['steps'] = None
             if reset_full:
                 _cache['full'] = None
                 _cache['hasdups'] = None
         return r

def iscomplexobj(x):
    """
    Check for a complex type or an array of complex numbers.

    The type of the input is checked, not the value. Even if the input
    has an imaginary part equal to zero, `iscomplexobj` evaluates to True.

    Parameters
    ----------
    x : any
        The input can be of any type and shape.

    Returns
    -------
    iscomplexobj : bool
        The return value, True if `x` is of a complex type or has at least
        one complex element.

    See Also
    --------
    isrealobj, iscomplex

    Examples
    --------
    >>> np.iscomplexobj(1)
    False
    >>> np.iscomplexobj(1+0j)
    True
    >>> np.iscomplexobj([3, 1+0j, True])
    True

    """
    return issubclass(asarray(x).dtype.type, _nx.complexfloating)

def guess_freq(dates):
    """Tries to estimate the frequency of a list of dates, by checking the steps
    between consecutive dates The steps should be in days.
    Returns a frequency code (alpha character)."""
    ddif = numeric.asarray(numpy.diff(dates))
    ddif.sort()
    if ddif.size == 0:
        fcode = _c.FR_UND
    elif ddif[0] == ddif[-1] == 1.:
        fcode = _c.FR_DAY
    elif (ddif[0] == 1.) and (ddif[-1] == 3.):
        fcode = _c.FR_BUS
    elif (ddif[0] > 3.) and  (ddif[-1] == 7.):
        fcode = _c.FR_WK
    elif (ddif[0] >= 28.) and (ddif[-1] <= 31.):
        fcode = _c.FR_MTH
    elif (ddif[0] >= 90.) and (ddif[-1] <= 92.):
        fcode = _c.FR_QTR
    elif (ddif[0] >= 365.) and (ddif[-1] <= 366.):
        fcode = _c.FR_ANN
    elif numpy.abs(24.*ddif[0] - 1) <= 1e-5 and \
         numpy.abs(24.*ddif[-1] - 1) <= 1e-5:
        fcode = _c.FR_HR
    elif numpy.abs(1440.*ddif[0] - 1) <= 1e-5 and \
         numpy.abs(1440.*ddif[-1] - 1) <= 1e-5:
        fcode = _c.FR_MIN
    elif numpy.abs(86400.*ddif[0] - 1) <= 1e-5 and \
         numpy.abs(86400.*ddif[-1] - 1) <= 1e-5:
        fcode = _c.FR_SEC
    else:
        warnings.warn("Unable to estimate the frequency! %.3f<>%.3f" %\
                      (ddif[0], ddif[-1]))
        fcode = _c.FR_UND
    return fcode