Python numeric.arange函数代码示例

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 testNoisyLine1(self):
   x = map(lambda x: x + gauss(0,0.002), arange(-1,1,0.001))
   y = map(lambda x: x + gauss(0,0.002), arange(-1,1,0.001))
   z = map(lambda x: x + gauss(0,0.02), arange(-1,1,0.001))
   line = array(zip(x,y,z))
   lpc = LPCImpl(h = 0.2, mult = 2)
   lpc_curve = lpc.lpc(X = line)

 def testNoisyLine2(self):
   x = map(lambda x: x + gauss(0,0.005), arange(-1,1,0.005))
   y = map(lambda x: x + gauss(0,0.005), arange(-1,1,0.005))
   z = map(lambda x: x + gauss(0,0.005), arange(-1,1,0.005))
   line = array(zip(x,y,z))
   lpc = LPCImpl(h = 0.2, convergence_at = 0.001, mult = 2)
   lpc_curve = lpc.lpc(X = line) 

def plot2():
  fig5 = plt.figure()
  x = map(lambda x: x + gauss(0,0.02)*(1-x*x), arange(-1,1,0.001))
  y = map(lambda x: x + gauss(0,0.02)*(1-x*x), arange(-1,1,0.001))
  z = map(lambda x: x + gauss(0,0.02)*(1-x*x), arange(-1,1,0.001))
  line = array(zip(x,y,z))
  lpc = LPCImpl(h = 0.05, mult = 2, it = 200, cross = False, scaled = False, convergence_at = 0.001)
  lpc_curve = lpc.lpc(X=line)
  ax = Axes3D(fig5)
  ax.set_title('testNoisyLine2')
  curve = lpc_curve[0]['save_xd']
  ax.scatter(x,y,z, c = 'red')
  ax.plot(curve[:,0],curve[:,1],curve[:,2])
  saveToPdf(fig5, '/tmp/testNoisyLine2.pdf')
  residuals_calc = LPCResiduals(line, tube_radius = 0.05, k = 10)
  residual_diags = residuals_calc.getPathResidualDiags(lpc_curve[0])
  fig6 = plt.figure()
  #plt.plot(lpc_curve[0]['lamb'][1:], residual_diags['line_seg_num_NN'], drawstyle = 'step', linestyle = '--')
  plt.plot(lpc_curve[0]['lamb'][1:], residual_diags['line_seg_mean_NN'])
  plt.plot(lpc_curve[0]['lamb'][1:], residual_diags['line_seg_std_NN'])
  saveToPdf(fig6, '/tmp/testNoisyLine2PathResiduals.pdf')
  coverage_graph = residuals_calc.getCoverageGraph(lpc_curve[0], arange(0.001, .102, 0.005))
  fig7 = plt.figure()
  plt.plot(coverage_graph[0],coverage_graph[1])
  saveToPdf(fig7, '/tmp/testNoisyLine2Coverage.pdf')
  residual_graph = residuals_calc.getGlobalResiduals(lpc_curve[0])
  fig8 = plt.figure()
  plt.plot(residual_graph[0], residual_graph[1])
  saveToPdf(fig8, '/tmp/testNoisyLine2Residuals.pdf')
  fig9 = plt.figure()
  plt.plot(range(len(lpc_curve[0]['lamb'])), lpc_curve[0]['lamb'])
  saveToPdf(fig9, '/tmp/testNoisyLine2PathLength.pdf')

def helixHeteroscedasticDiags():
  #Parameterise a helix (no noise)
  fig5 = plt.figure()
  t = arange(-1,1,0.0005)
  x = map(lambda x: x + gauss(0,0.001 + 0.001*sin(2*pi*x)**2), (1 - t*t)*sin(4*pi*t))
  y = map(lambda x: x + gauss(0,0.001 + 0.001*sin(2*pi*x)**2), (1 - t*t)*cos(4*pi*t))
  z = map(lambda x: x + gauss(0,0.001 + 0.001*sin(2*pi*x)**2), t)
  line = array(zip(x,y,z))
  lpc = LPCImpl(h = 0.1, t0 = 0.1, mult = 1, it = 500, scaled = False, cross = False)
  lpc_curve = lpc.lpc(X=line)
  ax = Axes3D(fig5)
  ax.set_title('helixHeteroscedastic')
  curve = lpc_curve[0]['save_xd']
  ax.scatter(x,y,z, c = 'red')
  ax.plot(curve[:,0],curve[:,1],curve[:,2])
  saveToPdf(fig5, '/tmp/helixHeteroscedastic.pdf')
  residuals_calc = LPCResiduals(line, tube_radius = 0.2, k = 20)
  residual_diags = residuals_calc.getPathResidualDiags(lpc_curve[0])
  fig6 = plt.figure()
  #plt.plot(lpc_curve[0]['lamb'][1:], residual_diags['line_seg_num_NN'], drawstyle = 'step', linestyle = '--')
  plt.plot(lpc_curve[0]['lamb'][1:], residual_diags['line_seg_mean_NN'])
  plt.plot(lpc_curve[0]['lamb'][1:], residual_diags['line_seg_std_NN'])
  saveToPdf(fig6, '/tmp/helixHeteroscedasticPathResiduals.pdf')
  coverage_graph = residuals_calc.getCoverageGraph(lpc_curve[0], arange(0.01, .052, 0.01))
  fig7 = plt.figure()
  plt.plot(coverage_graph[0],coverage_graph[1])
  saveToPdf(fig7, '/tmp/helixHeteroscedasticCoverage.pdf')
  residual_graph = residuals_calc.getGlobalResiduals(lpc_curve[0])
  fig8 = plt.figure()
  plt.plot(residual_graph[0], residual_graph[1])
  saveToPdf(fig8, '/tmp/helixHeteroscedasticResiduals.pdf')
  fig9 = plt.figure()
  plt.plot(range(len(lpc_curve[0]['lamb'])), lpc_curve[0]['lamb'])
  saveToPdf(fig9, '/tmp/helixHeteroscedasticPathLength.pdf')

    def GetFlow(self):
        '''
        Calculating inlet flow (coefficients of the FFT  x(t)=A0+sum(2*Ck*exp(j*k*2*pi*f*t)))
        Timestep and period from SimulationContext are necessary.
        '''
        try:
            timestep = self.SimulationContext.Context['timestep']
        except KeyError:
            print "Error, Please set timestep in Simulation Context XML File"
            raise
        try:
            period = self.SimulationContext.Context['period']
        except KeyError:
            print "Error, Please set period in Simulation Context XML File"
            raise

        t = arange(0.0,period+timestep,timestep).reshape((1,ceil(period/timestep+1.0)))
        Cc = self.f_coeff*1.0/2.0*1e-6
        Flow = zeros((1, ceil(period/timestep+1.0)))
        for freq in arange(0,ceil(period/timestep+1.0)):
            Flow[0, freq] = self.A0_v
            for k in arange(0,self.f_coeff.shape[0]):
                Flow[0, freq] = Flow[0, freq]+real(2.0*complex(Cc[k,0],Cc[k,1])*exp(1j*(k+1)*2.0*pi*t[0,freq]/period))
        self.Flow = Flow
        return Flow

    def GetTimeFlow(self, el, time):
        '''
        Calculating inlet flow (coefficients of the FFT  x(t)=A0+sum(2*Ck*exp(j*k*2*pi*f*t)))
        for a specific time value.
        If signal is specified, flow is computed from time values.
        '''
        try:
            period = self.SimulationContext.Context['period']
        except KeyError:
            print "Error, Please set period in Simulation Context XML File"
            raise

        try:
            signal = self.InFlows[el]['signal']
            try:
                timestep = self.SimulationContext.Context['timestep']
            except KeyError:
                print "Error, Please set timestep in Simulation Context XML File"
                raise
            t = arange(0.0,period+timestep,timestep)
            t2 = list(t)
            Flow = float(signal[t2.index(time)])/6.0e7
            self.Flow = Flow
            return Flow
        except KeyError:
            f_coeff = self.InFlows[el]['f_coeff']
            A0 = self.InFlows[el]['A0']
            Cc = f_coeff*1.0/2.0*1e-6
            Flow = A0
            for k in arange(0,f_coeff.shape[0]):
                Flow += real(2.0*complex(Cc[k,0],Cc[k,1])*exp(1j*(k+1)*2.0*pi*time/period))
            self.Flow = Flow
            return Flow

 def __getitem__(self, key):
     try:
         size = []
         typ = int
         for k in range(len(key)):
             step = key[k].step
             start = key[k].start
             if start is None:
                 start = 0
             if step is None:
                 step = 1
             if isinstance(step, complex):
                 size.append(int(abs(step)))
                 typ = float
             else:
                 size.append(
                     int(math.ceil((key[k].stop - start)/(step*1.0))))
             if (isinstance(step, float) or
                     isinstance(start, float) or
                     isinstance(key[k].stop, float)):
                 typ = float
         if self.sparse:
             nn = [_nx.arange(_x, dtype=_t)
                     for _x, _t in zip(size, (typ,)*len(size))]
         else:
             nn = _nx.indices(size, typ)
         for k in range(len(size)):
             step = key[k].step
             start = key[k].start
             if start is None:
                 start = 0
             if step is None:
                 step = 1
             if isinstance(step, complex):
                 step = int(abs(step))
                 if step != 1:
                     step = (key[k].stop - start)/float(step-1)
             nn[k] = (nn[k]*step+start)
         if self.sparse:
             slobj = [_nx.newaxis]*len(size)
             for k in range(len(size)):
                 slobj[k] = slice(None, None)
                 nn[k] = nn[k][slobj]
                 slobj[k] = _nx.newaxis
         return nn
     except (IndexError, TypeError):
         step = key.step
         stop = key.stop
         start = key.start
         if start is None:
             start = 0
         if isinstance(step, complex):
             step = abs(step)
             length = int(step)
             if step != 1:
                 step = (key.stop-start)/float(step-1)
             stop = key.stop + step
             return _nx.arange(0, length, 1, float)*step + start
         else:
             return _nx.arange(start, stop, step)

def tri(N, M=None, k=0, dtype=float):
    """ returns a N-by-M array where all the diagonals starting from
        lower left corner up to the k-th are all ones.
    """
    if M is None: M = N
    m = greater_equal(subtract.outer(arange(N), arange(M)),-k)
    return m.astype(dtype)

def eye(N, M=None, k=0, dtype=float):
    """ eye returns a N-by-M 2-d array where the  k-th diagonal is all ones,
        and everything else is zeros.
    """
    if M is None: M = N
    m = equal(subtract.outer(arange(N), arange(M)),-k)
    if m.dtype != dtype:
        return m.astype(dtype)

def diagflat(v, k=0):
    """
    Create a two-dimensional array with the flattened input as a diagonal.

    Parameters
    ----------
    v : array_like
        Input data, which is flattened and set as the `k`-th
        diagonal of the output.
    k : int, optional
        Diagonal to set; 0, the default, corresponds to the "main" diagonal,
        a positive (negative) `k` giving the number of the diagonal above
        (below) the main.

    Returns
    -------
    out : ndarray
        The 2-D output array.

    See Also
    --------
    diag : MATLAB work-alike for 1-D and 2-D arrays.
    diagonal : Return specified diagonals.
    trace : Sum along diagonals.

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

    >>> np.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):
    """
    Extract a diagonal or construct a diagonal array.

    Parameters
    ----------
    v : array_like
        If `v` is a 2-dimensional array, return a copy of
        its `k`-th diagonal. If `v` is a 1-dimensional array,
        return a 2-dimensional array with `v` on the `k`-th diagonal.
    k : int, optional
        Diagonal in question.  The defaults is 0.

    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(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 = 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 plot1():
  fig1 = plt.figure()
  x = map(lambda x: x + gauss(0,0.005), arange(-1,1,0.005))
  y = map(lambda x: x + gauss(0,0.005), arange(-1,1,0.005))
  z = map(lambda x: x + gauss(0,0.005), arange(-1,1,0.005))
  line = array(zip(x,y,z))
  lpc = LPCImpl(h = 0.05, mult = 2, scaled = False)
  lpc_curve = lpc.lpc(X=line)
  ax = Axes3D(fig1)
  ax.set_title('testNoisyLine1')
  curve = lpc_curve[0]['save_xd']
  ax.scatter(curve[:,0],curve[:,1],curve[:,2],c = 'red')
  return fig1

def tri(N, M=None, k=0, dtype=float):
    """
    An array with ones at and below the given diagonal and zeros elsewhere.

    Parameters
    ----------
    N : int
        Number of rows in the array.
    M : int, optional
        Number of columns in the array.
        By default, `M` is taken equal to `N`.
    k : int, optional
        The sub-diagonal at and below which the array is filled.
        `k` = 0 is the main diagonal, while `k` < 0 is below it,
        and `k` > 0 is above.  The default is 0.
    dtype : dtype, optional
        Data type of the returned array.  The default is float.


    Returns
    -------
    tri : ndarray of shape (N, M)
        Array with its lower triangle filled with ones and zero elsewhere;
        in other words ``T[i,j] == 1`` for ``i <= j + k``, 0 otherwise.

    Examples
    --------
    >>> np.tri(3, 5, 2, dtype=int)
    array([[1, 1, 1, 0, 0],
           [1, 1, 1, 1, 0],
           [1, 1, 1, 1, 1]])

    >>> np.tri(3, 5, -1)
    array([[ 0.,  0.,  0.,  0.,  0.],
           [ 1.,  0.,  0.,  0.,  0.],
           [ 1.,  1.,  0.,  0.,  0.]])

    """
    if M is None:
        M = N

    m = greater_equal.outer(arange(N, dtype=_min_int(0, N)),
                            arange(-k, M-k, dtype=_min_int(-k, M - k)))

    # Avoid making a copy if the requested type is already bool
    if np_dtype(dtype) != np_dtype(bool):
        m = m.astype(dtype)

    return m