def wiki_code():
    N = 51
    # Create a sparse matrix in the "list of lists" format. 
    # This provides a flexible syntax for creating sparse matrices.
    A = scpy.sparse.lil_matrix((N, N))
    # Finite difference matrices will always have some 
    # regular structure involving diagonals and off diagonals.
    # The lil_matrix object has a function to help you with this.
    A.setdiag(np.ones(N))              # The diagonal
    A.setdiag(-2*np.ones(N-1),k=1)       # The fist upward off-diagonal.
    A.setdiag(2*np.ones(N-1), k=-1)
    # etc. observe that I leave it to you to insert correct values along the diagonals...

    # To fix the boundaries, you'll have to reset some rows.
    # The numpy indexing you expect to work will work.
    #A[0,:] = np.zeros(N) # zero out the first row
    #A[0,0] = 1.0       # set diagonal in that row to 1.0


    # For performance, other sparse matrix formats are preferred.
    # Convert to "Compressed Row Format" CR
    A = A.tocsr()

    # Helpful for diagnostics

    print A.todense()

    # and 

    plt.spy(A.todense())
    plt.show()

def test_tsqr(create_func):
    mat, data = create_func()
    n = mat.shape[1]

    q, r = csnmf.tsqr.qr(data)

    print q.shape
    q = np.array(q)

    r = np.array(r)
    print r.shape

    print np.linalg.norm(mat - np.dot(q, r))

    assert np.allclose(mat, np.dot(q, r))
    assert np.allclose(np.eye(n, n), np.dot(q.T, q))
    assert np.all(r == np.triu(r))

    plt.figure()
    plt.subplot(2, 4, 1)
    plt.imshow(mat, interpolation='nearest')
    plt.title('Original matrix')
    plt.subplot(2, 4, 2)
    plt.imshow(q, interpolation='nearest')
    plt.title('$\mathbf{Q}$')
    plt.subplot(2, 4, 3)
    plt.imshow(np.dot(q.T, q), interpolation='nearest')
    plt.title('$\mathbf{Q}^T \mathbf{Q}$')
    plt.subplot(2, 4, 4)
    plt.imshow(r, interpolation='nearest')
    plt.title('$\mathbf{R}$')

    plt.subplot(2, 4, 8)
    plt.spy(r)
    plt.title('Nonzeros in $\mathbf{R}$')

def plot_matrix_fig(matrix, filename):
	fig = plt.figure()
	plt.spy(matrix)
	plt.title(filename)
	#plt.savefig('/home/igorpesic/Desktop/filename.png', dpi=600)
	plt.show()
	return

def sparsity(matrix, min_value=0):
    """
    https://redmine.epfl.ch/projects/python_cookbook/wiki/Matrix_sparsity_patterns
    """
    mat = matrix.matrix
    plt.spy(mat, precision=min_value, marker=',')
    plt.show()

def plotFeatMatr(toPlot,featureObjects,featureMat,saveDir,label,badIdx):
    nFeats = featureMat.shape[1]
    nnzPerFeature = toPlot.getnnz(0)
    # get the indices to sort this ish.
    # how many should we use?...
    # get the top N most common
    mostCommon = np.argsort(nnzPerFeature)[-nFeats//7:]
    # get their labels
    featLabels = [f.label() for f in featureObjects]
    # get a version imshow can handle
    matImage = toPlot.todense()
    # fix the aspect ratio
    aspectSkew = len(badIdx)/nFeats
    aspectStr = 1./aspectSkew
    # plot everything
    ax = plt.subplot(1,1,1)
    cax = plt.imshow(matImage,cmap=plt.cm.hot_r,aspect=aspectStr,
                     interpolation="nearest")
    plt.spy(toPlot,marker='s',markersize=1.0,color='b',
            aspect=aspectStr,precision='present')
    cbar = plt.colorbar(cax, ticks=[0, 1], orientation='vertical')
    # horizontal colorbar
    cbar.ax.set_yticklabels(['Min Feat', 'Max Feat'],
                            fontsize=g_label)
    ax.set_xticks(range(nFeats))
    ax.set_xticklabels(featLabels,rotation='vertical')
    plt.xlabel("Feature Number",fontsize=g_label)
    plt.ylabel("Individual",fontsize=g_label)
    return aspectStr

def plot(f, mesh=None, fig=None, **kwargs):
    """Plot functions/expression in 1D and spy matrices."""
    # Decide if function or expression
    try:
        V = f.function_space
        # Not really Scattered
        x = V.mesh.vertex_coordinates
        # Scattered in the same way
        y = f.vertex_values()
        isort = np.argsort(x[:, 0])
        x = x[isort]
        y = y[isort]
    # Expression or matrix
    except AttributeError:
        try:
            F = f.toarray()
            fig = plt.figure()
            plt.spy(F, precision=1e-10)
            return fig
        # Expression
        except AttributeError:
            # Evaluate at mesh vertices
            assert mesh is not None
            x = np.sort(mesh.vertex_coordinates)
            y = f.eval(x)

    if fig is None:
        fig = plt.figure()
    ax = fig.gca()
    ax.plot(x, y, **kwargs)

    return fig

def save_matrix_fig(matrix, path, filename):
	fig = plt.figure()
	plt.spy(matrix)
	#plt.title(filename)
	plt.savefig(path + filename + '.png', dpi=600)
	plt.close(fig)
	return

def makespyplot( matrix, name, imgtype=None ):
  if not scipy.sparse.isspmatrix( matrix ):
      matrix = matrix.toscipy()

  with plot.PyPlot( name, ndigits=0, imgtype=imgtype ) as plt:
    plt.spy( matrix, markersize=0.8, color='black')
    plt.title( name+', nnz = '+str(matrix.nnz) )

def plotSpy(M, eps, myTitle):
    plt.figure()
    frame1 = plt.gca()
    frame1.get_xaxis().set_ticks([])
    frame1.get_yaxis().set_ticks([])
    #  plt.title(myTitle+",eps="+str(eps))
    plt.spy(M, precision=eps, marker=".", markersize=3)
    plt.savefig(myTitle + "eps" + str(eps) + ".eps")

    def getError(self):
        if plotIt:
            plt.spy(self.getAve(self.M))
            plt.show()

        num = self.getAve(self.M) * self.getHere(self.M)
        err = np.linalg.norm((self.getThere(self.M)-num), np.inf)
        return err

def plot_sparse(mat, filename):
    plt.spy(mat, markersize=1)
    plt.tight_layout()
    #plt.rc('xtick',labelsize=16)
    #plt.rc('ytick',labelsize=16)
    plt.savefig.format = "pdf"
    plt.savefig(filename)
    plt.clf()

def plotMatrix(A):
    ''' Spy plot of matrix A.
    Arguments:
        A: (numpy.array, numpy.matrix, sparse.matrix)
    '''
    plt.spy(A, markersize=5)    
    plt.show()
    return

def plotSparsity(data, xmin=0, ymin=0, xmax=0, ymax=0):
    
    if xmax == 0:
        xmax = data.shape[0]
    if ymax == 0:
        ymax = data.shape[1]
        
    plt.spy(data[xmin:xmax, ymin:ymax])
    plt.show()

def makeCouplingAtomPlot(n, tag):
    import parseXYZ
    from mpl_toolkits.mplot3d import Axes3D

    mmdir = "/Volumes/s/matrices/matrixmarket/"
    mmfile = mmdir + tag + "_A.mm"
    xyzdir = "/Volumes/s/keceli/Dropbox/work/SEMO/xyz/"
    xyzfile = xyzdir + tag + ".xyz"
    mat = getAtomMatrix(mmfile)
    nnz = mat.nnz
    distW = np.array([0.0] * nnz)
    xyz = parseXYZ.getXYZ(xyzfile)
    matD = getDistMat(xyz, 5.6)
    makeNnzHist(mat, xyz)
    xyzW = np.asarray(xyz[1:4])
    myrow = np.asarray(mat.row)
    mycol = np.asarray(mat.col)
    listofAtoms1 = myrow[mycol == n]
    listofAtoms2 = mycol[myrow == n]
    listofAtoms = np.unique(np.concatenate((listofAtoms1, listofAtoms2)))
    print "focused atom number and its coordinates:", n + 1, xyzW[:, n]
    print "number of atoms interacting with this atom:", len(listofAtoms)
    print "list of interacting atoms:", listofAtoms + 1
    print "jmol command:", str(["C" + str(listofAtoms[i] + 1) for i in range(len(listofAtoms))]).replace(
        "[", ""
    ).replace("]", "").replace("'", "")
    print "coordinates of interacting atoms", xyzW[:, listofAtoms]
    fig = plt.figure()
    plt.spy(mat, markersize=1)
    fig = plt.figure()
    plt.spy(matD, markersize=1)
    fig = plt.figure()
    plt.plot(listofAtoms, linestyle="None", marker="s")
    fig = plt.figure()
    ax = fig.add_subplot(111, projection="3d")
    plt.plot(xyzW[0, :], xyzW[1, :], xyzW[2, :], linestyle="None", marker=".", color="black", markersize=2)
    plt.plot(
        xyzW[0, listofAtoms],
        xyzW[1, listofAtoms],
        xyzW[2, listofAtoms],
        linestyle="None",
        marker="o",
        color="red",
        markersize=3,
    )
    plt.plot([xyzW[0, n]], [xyzW[1, n]], [xyzW[2, n]], linestyle="None", marker="s", color="blue", markersize=4)
    plt.grid(which="major", axis="all")
    plt.gca().set_aspect("equal", adjustable="box")
    #    ax.auto_scale_xyz()
    # plt.axis('equal')
    ax.set_xlim([-5, 25])
    ax.set_ylim([-5, 25])
    ax.set_zlim([0, 45])
    ax.set_xlabel("x")
    ax.set_ylabel("y")
    ax.set_zlabel("z")
    return

def contact_to_distance(M):
    N = np.zeros(M.shape)
    N[M!=0] = 1/M[M!=0]
    s = sparse.csgraph.floyd_warshall(N)
    plt.figure()
    plt.spy(s, precision=0, markersize=0.00000000000005)
    plot_image = plt.imshow(s, vmin=0,vmax=np.percentile(s,99),interpolation='nearest')
    plot_image.set_cmap("jet")
    plt.show(block=BLOCK)
    print s
    return s

def plotSpy(M, eps):
    plt.figure()
    frame1 = plt.gca()
    frame1.get_xaxis().set_ticks([])
    frame1.get_xaxis().set_ticklabels([])
    frame1.get_yaxis().set_ticks([])
    frame1.get_yaxis().set_ticklabels([])
    plt.subplot(111)
    #  plt.title(myTitle+",eps="+str(eps))
    plt.spy(M, precision=eps, marker=".", markersize=1)
    plt.show

def plot_matrix(Mat,dir_name=None,title=None):
	"""Plots and saves a dense matrix"""
	if dir_name is None:
		dir_name='./'
	if title is None:
		title = 'Matrix'
	plt.figure(figsize=(10,10))
	plt.spy(np.abs(Mat))
	plt.grid()
	plt.title(title)
	plt.savefig(dir_name+title+'.png')

def main():
	n = 9
	arr = np.zeros((n,n))

	arr[:,0:3] = 1
	arr[n-1, :]  = 1
	arr[(4,7,1),(5,7,8)] = 1

	plt.spy(arr)
	savefig('arr_plot.png')
	plt.show()

 def plot_Phis_sparsity(self, Phis, fig=0):
     Phis = [phis[0].data.cpu().numpy() for phis in Phis]
     plt.figure(fig)
     plt.clf()
     for i, phi in enumerate(Phis):
         plt.subplot(1, len(Phis), i + 1)
         # plot first element of the batch
         plt.spy(phi, precision=0.001, marker='o', markersize=2)
         plt.xticks([])
         plt.yticks([])
         plt.title('k={}'.format(i))
     path = os.path.join(self.path, 'Phis.png')
     plt.savefig(path)

    def PlotSparsity(self, solver):
        """
        Call in with your T.Solver as argument.
        """
        solver['H'].evaluate()
        H = solver['H'].output()

        solver['dg'].evaluate()
        dg = solver['dg'].output()
        
        plt.figure(1)
        plt.spy(np.concatenate([H,dg.T],axis=1))
        plt.show()