Commits

cls  committed 6184834

graph.api.Graph, graph.generator.*

  • Participants
  • Parent commits da7c829

Comments (0)

Files changed (5)

File src/example/lamg_example.py

 import os
 import logging
 import datetime
+import numpy
 
 
 from lamg.lin import factory
         #%---------------------------------------------------------------------
         #% Create a graph Laplacian
         #%---------------------------------------------------------------------
-        G = graphs.Graphs.grid("fd", (400, 400), normalized=True)
+        G = graphs.Graphs.grid("fd", numpy.array([400, 400]), normalized=True)
         
         #%---------------------------------------------------------------------
         #% Setup phase: construct a LAMG multi-level hierarchy

File src/graph/api.py

 import networkx
 import numpy
 import scipy.sparse
+import logging
 
-import mat
+from mat import lab
 from aux.lazy import LazyProperty
 
 
         # TODO: indices consistent with adjacency/laplacian
         degreeArray = numpy.array([degreeMap[i] for i in range(n)])
         return degreeArray
-        
+        
+        
+
+
+
+class Graph:
+    """
+        %GRAPH A weighted directed graph (network) data structure.
+    %   This class represents an immutable weighted undirected graph
+    %   G=(V,E) with nodes V and edges E between them. Each edge (u,v) has
+    %   an associated positive weight C(u,v).
+    %
+    %   We assume that nodes are not self-connected, and that there are no
+    %   double edges (i.e. if (u,v) is an edge, then (v,u) is not).
+    %   Therefore this class models both such direct graphs and undirect
+    %   graphs with fixed edge orientations.
+    %
+    %   See also: GRAPHPLOTTER, GRAPHLOADER.
+    """
+    
+    #    %=========================== PROPERTIES ==============================
+    #    properties (GetAccess = public, SetAccess = public)
+    #        coord               % Optional d-D node locations for graph drawing
+    #    end
+    #    properties (Dependent)
+    #        numNodes            % # nodes
+    #        numEdges            % # edges
+    #        edge                % Edge list (u,v,edgeId). edgeId is a running ID from 1 to numEdges
+    #        edgesAndWeights     % Edge + weight list (u,v,w(u,v))
+    #%        weight              % Edge weight list (w)
+    #%        weightMatrix        % Edge weight matrix (diag(w))
+    #%        incidence           % Incidence matrix
+    #    end
+    #    properties (GetAccess = public, SetAccess = private)
+    #        metadata            % Data about this graph
+    #        
+    #        % Cached dependent properties
+    #%        adjacency           % Undirected adjacency matrix. Always upper-diagonal.
+    #        adjacency        % Symmetrized adjacency matrix (adjacency+adjacency')
+    #        degree              % Node degree list
+    #        laplacian           % Laplacian matrix
+    #        mass                % Laplacian diagonal mass matrix
+    #%       orientation         % Conventional edge orientation matrix
+    #%        edgeIndex           % Symmetrized matrix with adjacency matrix non-zero pattern. Entry(i,j) = row index of this edge in the obj.edge list
+    #    end
+    
+    
+    def __init__(self, metadata, dataType, data, coord):
+        """
+                    %Graph Constructor.
+            %   Graph(metadata, A, coord) constructs a graph with specified
+            %   metadata, symmetric adjacency matrix with no diagonal, and
+            %   optional d-D node locations for graph drawing.
+            %
+            %   If A is empty, it is ignored.
+        """
+        # % Read adjacency matrix
+        if (dataType is "sym-adjacency"):
+            # % Symmetric adjacency matrix
+            A = data
+        elif (dataType is "adjacency"):
+            # % data = upper-triangular part of the adjacency matrix.
+            # % Symmetrize and remove diagonal.
+            W = lab.triu(data, 1)
+            A = W + lab.ctranspose(W) #  % Seems to be slightly faster than max(W,W')
+        elif (dataType is "edge-list"):
+            # % data is a matrix whose row format is assumed to be [i
+            # % j aij].
+            maxData = lab.max_(data, None, 1)
+            numNodes = lab.max_(maxData[:, 0:1]) # CHECK: range bounds
+            W = lab.sparse(data[:, 0], data[:, 1], data[:3], numNodes, numNodes) # CHECK: indices
+            A = W + lab.ctranspose(W) #  % Seems to be slightly faster than max(W,W')
+        elif (dataType is "laplacian"):
+            # % data = graph Laplacian. Is NOT saved in the object.
+            A = data
+        else:
+            raise Exception("Unrecognized graph data type: %s" % dataType)
+        
+        # store metadata
+        self.metadata = metadata
+        metadata.graphType = graph.api.GraphType.UNDIRECTED #% Graph is always UNDIRECTED here
+        
+        if (dataType is "laplacian"):
+            #% data = graph Laplacian. Assumed to be singly-connected.
+            #% Is not saved in the object after computing and caching
+            #% dependent fields.
+            self.coord = coord
+            self.metadata.numNodes = lab.size(data, 1)
+            if (lab.size(data, 1) == 1):
+                self.metadata.numEdges = 0
+            else:
+                self.metadata.numEdges = (lab.nnz(data) - self.metadata.numNodes) / 2 # % Assuming zero diagonal!
+                
+            # % Cache matrix fields
+            self.laplacian = A
+            self.adjacency = lab.diag(lab.diag(A)) - A
+            self.degree = lab.sum_(data != 0, 1) - 1 # TODO: % Could be improved a little by MEX
+        else:
+            self.metadata.numNodes = lab.size(A, 1)
+            if (lab.size(A, 1) == 1):
+                self.metadata.numEdges = 0
+            else:
+                self.metadata.numEdges = lab.nnz(A) / 2 # % Assuming zero diagonal!
+                
+            # % Cache matrix fields
+            self.adjacency = A
+            self.degree = lab.sum_(A != 0, 1) # TODO: % Could be improved a little by MEX
+            
+        # % Read coordinates
+        if ((coord is not None) and (lab.size(coord, 1) != lab.size(A, 1))):
+            logging.getLogger("%s.Graph" % __name__).warning("Coordinate vector was not of size numNodes x d, ignoring")
+        else:
+            self.coord = coord
+        
+        
+    @staticmethod
+    def newNamedInstance(name, dataType, data, coord):
+        """
+            % newNamedInstance(name, numNodes, edgeData, coord)
+            % constructs a named graph with the specified name, edge matrix
+            % (edgeData(i,:)=(u,v,w) represents the
+            %   ith edge nodes u and v, carrying weight w) and optional 2-D
+            %   node locations for graph drawing.
+        """
+        md = GraphMetadata()
+        md.name = name
+        graph = Graph.newInstanceFromMetadata(md, dataType, data, coord)
+        return graph
+    
+    @staticmethod
+    def newInstanceFromMetadata(md, dataType, data, coord):
+        """
+            % newNamedInstance(name, numNodes, edgeData, coord)
+            % constructs a named graph with the specified name, edge matrix
+            % (edgeData(i,:)=(u,v,w) represents the
+            %   ith edge nodes u and v, carrying weight w) and optional 2-D
+            %   node locations for graph drawing.
+        """
+        graph = Graph(md, dataType, data, coord)
+        return graph
+    
+    
+    def subgraph(self, nodes, name=None):
+        """
+        % Return the subgraph of g containing the nodes "nodes" and
+            % their links to themselves only.
+        """
+        m = GraphMetadata.copy(self.metadata)
+        m.group = None # FIXME: [obj.metadata.group '/' obj.metadata.name]
+        if (name is None):
+            name = "subgraph"
+        m.name = name
+        m.file = None
+        if (self.coord is not None):
+            subCoord = self.coord[nodes, :]
+        else:
+            subCoord = None
+        G = Graph(m, "adjacency", self.adjacency[nodes, nodes])
+        return G
+    
+    
+    def edgeFunction(self, fdata, **kwargs):
+        """
+                    % Convert an array FDATA to an edge function F on the graph.
+            % FDATA must be of length NUMEDGES. F=OBJ.EDGEFUNCTION(FDATA)
+            % is a sparse matrix with the non-zero pattern of
+            % OBJ.ADJACENCY. F=OBJ.EDGEFUNCTION(FDATA,
+            % graph.api.GraphType.UNDIRECTED) is its symmetrized version
+            % (FDATA should still be of length NUMEDGES!).
+        """
+        raise NotImplementedError("TODO:")
+    
+    
+    def strongAdjacency(self, threshold):
+        """
+            % Return the filtered strong connection adjacency matrix for
+            % threshold THRESHOLD.
+        """
+        raise NotImplementedError("TODO:")
+    
+    
+    def strrongEdgePortion(self, threshold):
+        """
+             % Return the percentage of strong connetions in the adjacency
+            % matrix.
+        """
+        raise NotImplementedError("TODO:")
+    
+    
+    def weakEdgePortion(self, threshold):
+        """
+        % Return the percentage of weak (small) connetions in the
+            % adjacency matrix.
+        """
+        r = 1 - self.strrongEdgePortion(threshold)
+        return r
+    
+    
+    @property
+    def numNodes(self):
+        """
+        % Return the number of graph nodes.
+        """
+        return self.metadata.numNodes
+    
+    
+    @property
+    def numEdges(self):
+        """
+        % Return the number of graph edges.
+        """
+        return self.metadata.numEdges
+    
+    
+    @property
+    def edge(self):
+        """
+        % Returns the edge list.
+        """
+        raise NotImplementedError("TODO:")
+    
+    @property
+    def edgesAndWeights(self):
+        """
+         % Returns the edge list.
+        """
+        raise NotImplementedError("TODO:")
+    
+    
+    @property
+    def mass(self):
+        """
+        % Laplacian mass matrix
+        """
+        D = lab.spdiags(lab.sum_(self.adjacency, 2), 0, self.numNodes, self.numNodes)
+        return D
+    
+    
+    @property
+    def laplacian(self):
+        """
+        % Return the (sparse) graph Laplacian L = D-A. Note that L may
+            % not be SPD if the weights are negative. Cached.
+        """
+        if (self.laplacian is None):
+            A = self.adjacency
+            D = self.mass
+            self.laplacian = D - A
+        return self.laplacian
+        
+        
+
+            
+            
+            
+            
+            

File src/graph/generator.py

 '''
 
 from lamg import amg
+import graph
+from mat import lab
 
 
 class GeneratorFactory:
     
     @staticmethod
     def parseArgs(type, **kwargs):
-        raise NotImplementedError("TODO:")
-    
+        options = amg.api.Options()
+        
+        options.addOption("type", type, lambda x: x in ["union", "path", "grid", "sun"])
+        options.addOption("gridType", None, amg.api.Options.isString)
+        options.addOption("bc", "neumann", lambda x: x in ["neumann", "periodic"])
+        options.addOption("h", None, amg.api.Options.isNumeric)
+        options.addOption("n", None, amg.api.Options.isPositiveInteger, "Grid meshsize")
+        options.addOption("e", 1.0, amg.api.Options.isNumeric, "Anisotropy parameter")
+        options.addOption("normalized", False, amg.api.Options.isLogical, "Normalize Laplacian or not")
+        options.addOption("stencil", 1.0, amg.api.Options.isNumeric, "Grid stencil")
+        options.addOption("alpha", None, amg.api.Options.isNumeric, "rotation angle")
+        options.addOption("eps", None, lambda x: x >= 0, "anisotropy coefficient")
+        options.addOption("extraEdgeWeight", 1.0, amg.api.Options.isNumeric)
+        options.addOption("g1", None, amg.api.Options.isA(graph.api.Graph))
+        options.addOption("g2", None, amg.api.Options.isA(graph.api.Graph))
+        
+        overridden = amg.api.Options.fromKeywordArgs(default=options, **kwargs)
+        return overridden
+
+
     
     
 class AbstractGeneratorGrid(amg.api.Builder):
     %   See also: GRAPH.
     """
     
+#        dim         % Grid dimension
+#        n           % Grid size
+#        h           % Grid meshsize
+#        numNodes    % Total number of gridpoints
+#        edges       % Graph adjacency list
+#        normalized  % If true, normalizes the Laplacian to O(1) L1 row sums
+#        type        % Grid discretization type
+#        stencil     % Laplacian stencil. Row format: [stencil_nbhr_relative_index(1:d) coef]
+#        nbhrFinder  % Sets boundary conditions
+#        NBHR_FINDER_FACTORY = graph.generator.NbhrFinderFactory
+
+    
     def __init__(self, options, type):
         """
         """
-        raise NotImplementedError("TODO:")
-    
+        if (options.n is None):
+            raise Exception("Must specify grid size")
+        
+        self.n = options.n
+        self.normalized = options.normalized
+        
+        if (options.h is None):
+            # % No meshsize specified, use h = 1/n in all directions so
+            # % that domain = unit cube
+            self.h = lab.rdivide(1.0, options.n)
+        else:
+            # % Meshsize specified
+            self.h = options.h
+        
+        # % Initializations
+        if (lab.size(self.h, 0) > 1):
+            options.h = lab.ctranspose(options.h)
+            self.h = lab.ctranspose(options.h)
+            
+        if (lab.size(self.n, 0) > 1):
+            options.n = lab.ctranspose(options.n)
+            self.n = lab.ctranspose(self.n)
+            
+        self.dim = lab.numel(self.n)
+        self.numNodes = lab.prod(self.n)
+        
+        self.type = type
+        self.stencil = self.buildStencil(options)
+        
+        options.dim = self.dim
+        self.nbhrFinder = graph.generator.AbstractGeneratorGrid.NBHR_FINDER_FACTORY.newInstance(options)
+        
+        
+        
     
     def buildStencil(self, options):
         """
         % Build the Laplacian stencil from input options.
         """
-        raise NotImplementedError("TODO:")
+        raise NotImplementedError("abstract method")
     
     
     def buildCoord(self):

File src/lamg/amg/api.py

     @classmethod
     def hasType(cls, typ):
         return (lambda x: type(x) is typ)
+    isA = hasType
     
 
     def __init__(self):

File src/mat/lab.py

     of rows. If DIM > NDIMS(X), M will be 1.
     """
     # INFO: since python is 0-based, size(x,0) should return the number of rows
-    try:
+    if (type(A) is numpy.ndarray):
         if dim is None:
             return A.shape
         else:
             return A.shape[dim]
-    except:
+    elif (type(A) is PETSc.Mat):
         # petsc mat
         (m,n) = A.getSizes()[0]
         if dim is not None:
 
 
 
+def prod(X, dim=None):
+    """
+     prod Product of elements.
+    For vectors, prod(X) is the product of the elements of X. For
+    matrices, prod(X) is a row vector with the product over each
+    column. For N-D arrays, prod(X) operates on the first
+    non-singleton dimension.
+ 
+    prod(X,DIM) works along the dimension DIM. 
+ 
+    Example: If X = [0 1 2
+                     3 4 5]
+ 
+    then prod(X,1) is [0 4 10] and prod(X,2) is [ 0
+                                                 60]
+    """
+    if isvector(X):
+        return numpy.prod(X)
+    elif ismatrix(X):
+        raise NotImplementedError("TODO:")
+    else:
+        raise NotImplementedError()
+
+
+