def __init__(self, *shape):
     if len(shape) == 1 and isinstance(shape[0], tuple):
         shape = shape[0]
     x = as_strided(_nx.zeros(1), shape=shape,
                    strides=_nx.zeros_like(shape))
     self._it = _nx.nditer(x, flags=['multi_index', 'zerosize_ok'],
                           order='C')

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 _set_selection(self, val):
        oldval = self._selection
        self._selection = val

        datasource = getattr(self.plot, self.axis, None)

        if datasource is not None:

            mdname = self.metadata_name

            # Set the selection range on the datasource
            datasource.metadata[mdname] = val
            datasource.metadata_changed = {mdname: val}

            # Set the selection mask on the datasource
            selection_masks = \
                datasource.metadata.setdefault(self.mask_metadata_name, [])
            for index in range(len(selection_masks)):
#                if id(selection_masks[index]) == id(self._selection_mask):
                if True:
                    del selection_masks[index]
                    break

            # Set the selection mode on the datasource
            datasource.metadata[self.selection_mode_metadata_name] = \
                      self.selection_mode

            if val is not None:
                low, high = val
                data_pts = datasource.get_data()
                new_mask = (data_pts >= low) & (data_pts <= high)
                selection_masks.append(new_mask)
                self._selection_mask = new_mask
            else:
                # Set the selection mask to false.
                data_pts = datasource.get_data()
                new_mask = zeros_like(data_pts,dtype=bool)
                selection_masks.append(new_mask)
                self._selection_mask = new_mask

                
            datasource.metadata_changed = {self.mask_metadata_name: val}
        

        self.trait_property_changed("selection", oldval, val)

        for l in self.listeners:
            if hasattr(l, "set_value_selection"):
                l.set_value_selection(val)

        return

def fix(x, y=None):
    """
    Round to nearest integer towards zero.

    Round an array of floats element-wise to nearest integer towards zero.
    The rounded values are returned as floats.

    Parameters
    ----------
    x : array_like
        An array of floats to be rounded
    y : ndarray, optional
        Output array

    Returns
    -------
    out : ndarray of floats
        The array of rounded numbers

    See Also
    --------
    trunc, floor, ceil
    around : Round to given number of decimals

    Examples
    --------
    >>> np.fix(3.14)
    3.0
    >>> np.fix(3)
    3.0
    >>> np.fix([2.1, 2.9, -2.1, -2.9])
    array([ 2.,  2., -2., -2.])

    """
    x = nx.asanyarray(x)
    if y is None:
        y = nx.zeros_like(x)
    y1 = nx.floor(x)
    y2 = nx.ceil(x)
    y[...] = nx.where(x >= 0, y1, y2)
    return y

	def recommend(self, evidences, n):
		
		#Array of item bought
		evidencesVec = zeros_like(self.ItemPrior)
		evidencesVec[evidences] = 1
		
		p_z_newUser = self.pLSAmodel.folding_in(evidencesVec)
		
		p_item_z = self.pLSAmodel.p_w_z
		
		p_item_newUser = zeros((p_item_z.shape[0]))
		for i in range(p_item_z.shape[1]):
			p_item_newUser+=p_item_z[:,i] * p_z_newUser[0]
		
		
		ind = argsort(p_item_newUser)
		
		if n == -1:
			return ind
		else:
			return ind[-n:]		

def polyval(p, x):
    """Evaluate the polynomial p at x.  If x is a polynomial then composition.

    Description:

      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]

      x can be a sequence and p(x) will be returned for all elements of x.
      or x can be another polynomial and the composite polynomial p(x) will be
      returned.

      Notice:  This can produce inaccurate results for polynomials with
      significant variability. Use carefully.
    """
    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 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)` is returned for each element of `x`.
    If `x` is another polynomial then the composite polynomial `p(x(t))`
    is returned.

    Parameters
    ----------
    p : array_like or poly1d object
       1D array of polynomial coefficients (including coefficients equal
       to zero) from highest degree to the constant term, or an
       instance of poly1d.
    x : array_like or poly1d object
       A number, a 1D array of numbers, or an instance of poly1d, "at"
       which to evaluate `p`.

    Returns
    -------
    values : ndarray or poly1d
       If `x` is a poly1d instance, the result is the composition of the two
       polynomials, i.e., `x` is "substituted" in `p` and the simplified
       result is returned. In addition, the type of `x` - array_like or
       poly1d - governs the type of the output: `x` array_like => `values`
       array_like, `x` a poly1d object => `values` is also.

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

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

    References
    ----------
    .. [1] I. N. Bronshtein, K. A. Semendyayev, and K. A. Hirsch (Eng.
       trans. Ed.), *Handbook of Mathematics*, New York, Van Nostrand
       Reinhold Co., 1985, pg. 720.

    Examples
    --------
    >>> np.polyval([3,0,1], 5)  # 3 * 5**2 + 0 * 5**1 + 1
    76
    >>> np.polyval([3,0,1], np.poly1d(5))
    poly1d([ 76.])
    >>> np.polyval(np.poly1d([3,0,1]), 5)
    76
    >>> np.polyval(np.poly1d([3,0,1]), np.poly1d(5))
    poly1d([ 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 __init__(self, *shape):
     x = as_strided(_nx.zeros(1), shape=shape, strides=_nx.zeros_like(shape))
     self._it = _nx.nditer(x, flags=['multi_index'], order='C')

def main():
    Omega = 1
    bz = BelousovZhabotinskii(Omega)
    
    #u0 = bz.max_u()
    u0 = bz.u_stationary()*2
    v0 = bz.v_stationary()
    
    y0 = [u0,v0]
    
    max_t = 10.0
    t = np.arange(0,max_t,0.0001)
    u_nullcline = np.logspace(-4, 0, 100000)*Omega
    v_nullclines = bz.nullcline(u_nullcline)
    
    y = odeint(bz.dy_dt,y0,t)
    
    plt.Figure()
    plt.plot(u_nullcline, v_nullclines[0])
    plt.plot(u_nullcline, v_nullclines[1])
    plt.plot(y[:,0],y[:,1])
    plt.loglog()
    plt.xlim([5e-5*Omega,2e0*Omega])
    #plt.show()
    plt.savefig("bz_wave_phase_plot.png")
    
    plt.clf()
    plt.plot(t,y[:,0])
    plt.plot(t,y[:,1])
    plt.plot(t,bz.w(y[:,0],y[:,1]))
    plt.yscale('log')
    plt.xlim((0,5))
    #plt.show()
    plt.savefig("bz_wave_concen_versus_time.png")
    
    plt.clf()
    h = 0.001
    x = np.arange(0,20,h)
    u0 = zeros_like(x) + bz.u_stationary()
    v0 = zeros_like(x) + bz.v_stationary()
    u0[x<1] = bz.u_stationary()*2
    y = np.vstack((u0,v0))
    
    if 1:
        dt = 0.0000001
        iterations = int(max_t/dt)
        out_every = iterations/1000
        #out_every = 1000
        #plt.ion()
        plot_u, = plt.plot(x,u0)
        plot_v, = plt.plot(x,v0)
        plt.yscale('log')
        plt.ylim((bz.u_stationary()/10,bz.max_u()))
        #plt.show()
        dydt_old = bz.dy_dt_diffuse(y, t, h)
        for i in range(0,iterations):
            t = i*dt
        
            if (i%out_every == 0):
                plot_u.set_ydata(y[0,:])
                plot_v.set_ydata(y[1,:])
                #plt.draw()
                plt.savefig("bz_wave_" + '%04d'%i + ".png")
            
            dydt = bz.dy_dt_diffuse(y, t, h)
            #y = y + dt*dydt
            y = y + 3.0/2.0*dt*dydt - 0.5*dt*dydt_old
            dydt_old = dydt

 def laplace(self,u,h):
     ret = zeros_like(u)
     ret[1:-1] = (u[0:-2] - 2.0*u[1:-1] + u[2:])/h**2
     ret[0] = (-u[0] + u[1])/h**2
     ret[-1] = (u[-2] - u[-1])/h**2
     return ret