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

def plot_price_volume(buy_prices, buy_volumes, sell_prices, sell_volumes, show_intersections=False):
    
    intersect_x = None
    intersect_y = None
    
    for i in range(0, len(buy_prices)):
        
        x_buy = [r for r in buy_volumes[i]]
        y_buy = [r for r in buy_prices[i]]
        x_sell = [r for r in sell_volumes[i]]
        y_sell = [r for r in sell_prices[i]]
        f_buy = lambda x: interp(x, x_buy, y_buy)
        f_sell = lambda x: interp(x, x_sell, y_sell)
        
        if show_intersections == True:
   
            intersect_x, intersect_y = get_intersection_point(x_buy, y_buy, x_sell, y_sell, time=str(times[i]))
            plt.plot(intersect_x, intersect_y, 'ko')
            
        plt.plot(x_buy, f_buy(x_buy), 'k-', label='Demand Curve')#, linewidth=0.2, markersize=0.1)
        plt.plot(x_sell, f_sell(x_sell), 'k--', label='Supply Curve')
        #plt.annotate(str(intersect_x[0])+ ' | ' + str(intersect_y[0]), xy=(intersect_x, intersect_y),
        #             xytext=(intersect_x-3000, intersect_y-100))
        #plt.plot(41500, 80.16, 'ro')
    plt.xlabel("MWh")
    plt.ylabel("EUR/MWh")
    plt.legend(loc='best')
    plt.ylim(-200, 2000)
    plt.yticks(arange(-200, 2000, 200.0))
    plt.grid(axis='y')
    plt.autoscale()
    plt.show()

def create_dataset():
    dataset = SupervisedDataSet(1, 1)

    for x in arange(0, 4*pi, pi/30):
        dataset.addSample(x, sin(x))

    return dataset

 def checkLinearConsistence(self):
     '''
     This method checks and fixes the correct proportion between the meshes of each edge.
     '''
     for edge in self.GraphEdgeToMesh.iterkeys():
         if edge.Side =='venous':
             meshes = self.GraphEdgeToMesh[edge]      
             startingMesh = meshes[0]
             startingRadius = self.ElementIdsToElements[str(startingMesh.keys()[0])].Radius[0]
             endingMesh = meshes[len(meshes)-1]
             endingRadius = self.ElementIdsToElements[str(endingMesh.keys()[0])].Radius[len(self.ElementIdsToElements[str(endingMesh.keys()[0])].Radius)-1]
             elLen = edge.Length['value']/len(meshes)
             self.dz = elLen/1.0e5
             z = arange(0.0,elLen,self.dz)
             dr = (endingRadius-startingRadius)/(len(meshes))
             for mesh in meshes:
                 if mesh == startingMesh:
                     r1 = startingRadius
                     r2 = r1+dr
                 elif mesh == endingMesh:     
                     r1 = r2
                     r2 = endingRadius
                 else: 
                     r1 = r2
                     r2+=dr
                 r_z = r1+((dr/elLen)*z)
                 self.ElementIdsToElements[str(mesh.keys()[0])].Radius = r_z

def twoDisjointLinesWithMSClustering():
 
  t = arange(-1,1,0.002)
  x = map(lambda x: x + gauss(0,0.02)*(1-x*x), t)
  y = map(lambda x: x + gauss(0,0.02)*(1-x*x), t)
  z = map(lambda x: x + gauss(0,0.02)*(1-x*x), t)
  line1 = array(zip(x,y,z))
  line = vstack((line1, line1 + 3))
  lpc = LPCImpl(start_points_generator = lpcMeanShift(ms_h = 1), h = 0.05, mult = None, it = 200, cross = False, scaled = False, convergence_at = 0.001)
  lpc_curve = lpc.lpc(X=line)
  #Plot results
  fig = plt.figure()
  ax = Axes3D(fig)
  labels = lpc._startPointsGenerator._meanShift.labels_
  labels_unique = unique(labels)
  cluster_centers = lpc._startPointsGenerator._meanShift.cluster_centers_
  n_clusters = len(labels_unique)
  colors = cycle('bgrcmyk')
  for k, col in zip(range(n_clusters), colors):
    cluster_members = labels == k
    cluster_center = cluster_centers[k]
    ax.scatter(line[cluster_members, 0], line[cluster_members, 1], line[cluster_members, 2], c = col, alpha = 0.1)
    ax.scatter([cluster_center[0]], [cluster_center[1]], [cluster_center[2]], c = 'b', marker= '^')
    curve = lpc_curve[k]['save_xd']
    ax.plot(curve[:,0],curve[:,1],curve[:,2], c = col, linewidth = 3)
  plt.show()

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