Dr. David Harrington

6083 Reputation

21 Badges

19 years, 323 days
University of Victoria
Professor or university staff
Victoria, British Columbia, Canada

Social Networks and Content at Maplesoft.com

Maple Application Center
I am a professor of chemistry at the University of Victoria, BC, Canada, where my research areas are electrochemistry and surface science. I have been a user of Maple since about 1990.

MaplePrimes Activity

These are Posts that have been published by dharr

Recently Mapleprimes user @vs140580  asked here about finding the shortest returning walk in a graph (explained more below). I provided an answer using the labelled adjacency matrix. @Carl Love pointed out that the storage requirements are significant. They can be reduced by storing only the vertices in the walks and not the walks themselves. This can be done surprisingly easily by redefining how matrix multiplication is done using Maple's LinearAlgebra:-Generic package.

(The result is independent of the chosen vertex.)



Consider the following graph. We want to find, for a given vertex v, the shortest walk that visits all vertices and returns to v. A walk traverses the graph along the edges, and repeating an edge is allowed (as distinct from a path, in which an edge can be used only once).

G := AddVertex(PathGraph(5), [6, 7]); AddEdge(G, {{3, 7}, {4, 6}, {6, 7}}); DrawGraph(G, layout = circle, size = [250, 250])

GRAPHLN(undirected, unweighted, [1, 2, 3, 4, 5, 6, 7], Array(1..7, {(1) = {2}, (2) = {1, 3}, (3) = {2, 4}, (4) = {3, 5}, (5) = {4}, (6) = {}, (7) = {}}), `GRAPHLN/table/2`, 0)

n := NumberOfVertices(G)


A := AdjacencyMatrix(G)

As is well known, the (i, j) entry of A^k gives the number of walks of length k between vertices i and j. The labelled adjacency matrix shows the actual walks.

Alab := `~`[`*`](A, Matrix(n, n, symbol = omega))

Matrix(%id = 36893490455680637396)

For example, the (3,3) entry of Alab^2 has 3 terms, one for each walk of length 2 that starts and ends at vertex 3. The `ω__i,j` factors in a term give the edges visited.

So the algorithm to find the shortest walk is to raise Alab to successively higher powers until one of the terms in the diagonal entry for the vertex of interest has indices for all the vertices.


Matrix(%id = 36893490455709504684)

The expressions for higher powers get very large quickly, so an algorithm that only retains the sets of vertices in each term as a set of sets will use less storage space. So we can consider the following modified labelled adjacency matrix.

B := Matrix(n, n, proc (i, j) options operator, arrow; if A[i, j] = 1 then {{i, j}} else {} end if end proc)

Matrix(%id = 36893490455703852204)

Now we need to modify how we do matrix multiplication, but Maple has the LinearAlgebra:-Generic package to do this. We can redefine addition and multiplication to operate correctly on the sets of sets.

Consider two sets of sets of vertices for walks.

a := {{7}, {1, 3}, {2, 6, 7}}; b := {{1, 2}, {2, 3, 8}}

{{7}, {1, 3}, {2, 6, 7}}

{{1, 2}, {2, 3, 8}}

Addition is just combining the possibilities, and set union will do this. And addition of "zero" should add no possibilities, so we take {} as zero.

`union`(a, b); `union`(a, {})

{{7}, {1, 2}, {1, 3}, {2, 3, 8}, {2, 6, 7}}

{{7}, {1, 3}, {2, 6, 7}}

Multiplication is just combining all pairs by union. Notice here that {7} union {1,3,5} and {1,5} union {3,7} give the same result, but that we do not get duplicates in the set.

{seq(seq(`union`(i, j), `in`(i, a)), `in`(j, b))}

{{1, 2, 3}, {1, 2, 7}, {1, 2, 3, 8}, {1, 2, 6, 7}, {2, 3, 7, 8}, {2, 3, 6, 7, 8}}

The unit for multiplication should leave the set of sets unchanged, so {{}} can be used

{seq(seq(`union`(i, j), `in`(i, a)), `in`(j, {{}}))}

{{7}, {1, 3}, {2, 6, 7}}

And we should check that zero multiplied by a is zero

{seq(seq(`union`(i, j), `in`(i, a)), `in`(j, {}))}


Define these operations for the LinearAlgebraGeneric package:

F[`=`] := `=`; F[`/`] := `/`; F[`0`] := {}; F[`1`] := {{}}; F[`+`] := `union`; F[`*`] := proc (x, y) options operator, arrow; {seq(seq(`union`(i, j), `in`(i, x)), `in`(j, y))} end proc

Warning, (in F[*]) `j` is implicitly declared local

Warning, (in F[*]) `i` is implicitly declared local

Compare B^2 with Alab^2. We have lost information about the details of the walks except for the vertices visited.

LinearAlgebra:-Generic:-MatrixMatrixMultiply[F](B, B)

Matrix(%id = 36893490455680647164)

So here is a procedure for finding the length of the shortest walk starting and ending at vertex v.

  uses GraphTheory;
  local x,y,i,j,vv,A,B,F,n,vertset;
  if IsDirected(G) then error "graph must be undirected" end if;
  if not member(v,Vertices(G),'vv') then error "vertex not in graph" end if;
  if not IsConnected(G) then return infinity end if;
  F[`=`]:=`=`:F[`/`]:=`/`: # not required but have to be specified
  F[`*`]:=(x,y)->{seq(seq(i union j, i in x), j in y)};
  A:=Matrix(n,n, (i, j)-> if AdjacencyMatrix(G)[i, j] = 1 then {{i, j}} else {} end if);
  for i from 2 do
  until vertset in B[vv,vv];
end proc:

WalkLength(G, 1)



Download WalksGenericSetOfSets.mw

Recently here Mapleprimes member @Icz asked about generating all minimal edge cuts of a graph. I gave a brute force algorithm based on testng all bipartitions of the vertices. I then looked into improving the method by using fundamental cutsets to generate the cutsets. A description of the method and the ideas behind it is given here, together with a comparison of the two methods.

[Edit - see below for a newer version.]

Find all minimal edge cuts of a connected undirected graph using fundamental cutsets. (Assume no self-loops or multiple edges.)
The MinimalEdgeCuts procedure is in the startup code edit region (See Edit -> Startup Code.) (and is reproduced at the end of the worksheet).



Choose a graph.

G := Graph({{1, 2}, {1, 3}, {1, 4}, {2, 3}, {3, 4}}); DrawGraph(G, size = [200, 200])

GRAPHLN(undirected, unweighted, [1, 2, 3, 4], Array(1..4, {(1) = {2, 3, 4}, (2) = {1, 3}, (3) = {1, 2, 4}, (4) = {1, 3}}), `GRAPHLN/table/1`, 0)

The objective is to find all minimal edge cuts. If we partition the vertices of a graph in two non-empty sets (parts), then a cut is defined as a set of edges of the graph that have one end in the first part and the other end in the second part. Removal of these edges breaks the graph into two (or more) disconnected components.

For example the cut with edges {1,2} and {2,3} leaves two components with vertices [1,3,4] in one component and vertex [2] in the other component.

cut := {{1, 2}, {2, 3}}; IsCutSet(G, cut); DrawGraph(DeleteEdge(G, cut, inplace = false), size = [200, 200])

{{1, 2}, {2, 3}}


A cut is minimal if adding back any of the edges makes it connected again, i.e, the cut leave only two components. The above cut is minimal. Sometimes the term cutset is used to mean a minimal cut, but sometimes, as in Maple's IsCutSet, cutset means any removal of edges that increases the number of components.

This graph has 6 minimal cuts:

mincuts := MinimalEdgeCuts(G)

Vector[column](%id = 36893490386482949164)

This partitioning process suggests a brute force algorithm to find all minimal cuts. Find all partitions of the vertices into two (non-empty) sets. The two induced subgraphs G__1 and G__2 have some of the edges of the graph. All edges of the graph that are not in the two subgraphs form a cut. For example for the above partition we have:

partition := [[1, 3, 4], [2]]; G__1 := InducedSubgraph(G, partition[1]); 'Edges(G__1)' = Edges(G__1); G__2 := InducedSubgraph(G, partition[2]); 'Edges(G__2)' = Edges(G__2); cut := `minus`(`minus`(Edges(G), Edges(G__1)), Edges(G__2))

[[1, 3, 4], [2]]

GraphTheory:-Edges(G__1) = {{1, 3}, {1, 4}, {3, 4}}

GraphTheory:-Edges(G__2) = {}

{{1, 2}, {2, 3}}

Now we have to check that each component is connected, or that there are a total of two components. Since there are in this case, the cut is minimal.

IsConnected(G__1) and IsConnected(G__2); nops(ConnectedComponents(DeleteEdge(G, cut, inplace = false)))



Note that we do have to test that there are only two components. For example for the path graph 1--2--3, the partition [[2],[1,3]] has three components (the three individual vertices). The associated cut {{1,2},{2,3}} is not minimal because we could put edge {1,2} back and the graph would still be disconnected.


The number of partitions to be tested is 2^(n-1)-1, where n is the number of vertices. Any subset of the vertices could be in the first set (except the empty set and the set of all vertices), and then the other vertices must be in the second set. However that will give double the possibilities, since swapping set 1 and set 2 makes no difference. In the present case there are 2^3-1 = 7 possibilities, as seen below, where 0 indicates that the vertex in that position is in set 1 and 1 indicates that it is in set 2. For example, possibility 4 is the above partition [[1, 3, 4], [2]].

lprint("   1 2 3 4"); Print(Iterator:-SetPartitions(4, parts = 2), showrank)

   1 2 3 4
1: 0 0 0 1
2: 0 0 1 0
3: 0 0 1 1
4: 0 1 0 0
5: 0 1 0 1
6: 0 1 1 0
7: 0 1 1 1

The algorithm described is implemented as procedure MinimalEdgeCutsPartition in the startup code. (The minimal cuts are not produced in the same order as above.)


However, a better approach is by building up the cuts fron a set of fundamental cuts (or cutsets, using the word in the restrictive sense). These are associated with a spanning tree. We choose any spanning tree, which we highlight with red edges. This is a tree graph (connected with no cycles) that has the same vertices as the graph.

treeG := SpanningTree(G); HighlightEdges(G, treeG); DrawGraph(G, size = [200, 200])

GRAPHLN(undirected, unweighted, [1, 2, 3, 4], Array(1..4, {(1) = {2, 3, 4}, (2) = {1}, (3) = {1}, (4) = {1}}), `GRAPHLN/table/23`, 0)

Any spanning tree graph has n-1 edges. This is also the number of vertices minus the number of connected components (=1), which is the rank of the graph, which we will call r.
We will find one fundamental cutset for each edge of the tree, so there will be r fundamental cutsets.

NumberOfVertices(G)-1, NumberOfEdges(treeG), GraphRank(G)

3, 3, 3

Select one of the tree edges, say {1,2}. Removal of this edge partitions the tree into two components. We find the vertex partition and the corresponding cut by the procedure above. (We do not have to test for only two components because that follows from the tree structure.). This is the example we did above. The first three cuts in mincuts are the fundamental ones. Each one has one of the tree edges.

f[1], f[2], f[3] := entries(mincuts[1 .. 3], nolist)

{{1, 2}, {2, 3}}, {{1, 3}, {2, 3}, {3, 4}}, {{1, 4}, {3, 4}}

It is possible to generate all the cuts of the graph involving the tree edges by taking all 2^r = 8 subsets of the set of fundamental cutsets (this includes the empty set), and adding the elements of each subset with a ring sum operation.

A ring sum of two edge sets includes all edges in both sets except those that are common to both (the symmetric difference or disjunctive union):

`⊕` := proc (edges__1, edges__2) options operator, arrow; `minus`(`union`(edges__1, edges__2), `intersect`(edges__1, edges__2)) end proc

proc (edges__1, edges__2) options operator, arrow; `minus`(`union`(edges__1, edges__2), `intersect`(edges__1, edges__2)) end proc

(Maple has this as the symmdiff command, but not in infix form.)
We can also represent a set of edges by a vector, with entry 1 if the corresponding edge is in the set and 0 otherwise. The edge order we choose here is the order in Edges(G). Then the ring sum operation is the same as addition modulo 2 of the corresponding vectors, though we won't make much use of this in implementing the algorithm. Set up to convert to vectors.

edges := Edges(G); e := NumberOfEdges(G); edgetable := table(`~`[`=`](edges, [`$`(1 .. e)])); tovec := proc (edges) options operator, arrow; Vector(e, `~`[`=`](map(proc (x) options operator, arrow; edgetable[x] end proc, edges), 1)) end proc

{{1, 2}, {1, 3}, {1, 4}, {2, 3}, {3, 4}}

Let's try the ring sum on f[1] and f[2]. It is a cut, and since it now has two tree edges it is not a fundamental cut. (It happens to be minimal because it splits the graph into two components.)

f[1]*`⊕`*f[2] = `⊕`(f[1], f[2]); IsCutSet(G, rhs(%))

{{1, 2}, {2, 3}}*`⊕`*{{1, 3}, {2, 3}, {3, 4}} = {{1, 2}, {1, 3}, {3, 4}}


Do the same thing in the vector representation. The same logic about the tree edges (first three entries) shows the fundamental vectors are basis vectors and the three vectors V__1, V__2, `mod`(V__1+V__2, 2)are linearly independent (with respect to the addition mod 2 operation).

V[1], V[2], V[3] := tovec(f[1]), tovec(f[2]), tovec(f[3]); '`mod`(V[1]+V[2], 2)' = `mod`(V[1]+V[2], 2)

V[1], V[2], V[3] := Vector(5, {(1) = 1, (2) = 0, (3) = 0, (4) = 1, (5) = 0}), Vector(5, {(1) = 0, (2) = 1, (3) = 0, (4) = 1, (5) = 1}), Vector(5, {(1) = 0, (2) = 0, (3) = 1, (4) = 0, (5) = 1})

`mod`(V[1]+V[2], 2) = Vector[column](%id = 36893490386597633668)

The ring sum of any two cuts is another cut, which may be minimal (a cutset) or a disjoint union of cutsets (= two cutsets together). Combining the elements of the 8 subsets gives the eight cuts for this graph associated with the chosen tree graph, i.e., it has all combinations of tree edges. In vector form they are:

cutvecs := seq(seq(seq(`mod`(i*V[1]+j*V[2]+k*V[3], 2), `in`(i, [0, 1])), `in`(j, [0, 1])), `in`(k, [0, 1]))

Vector[column](%id = 36893490386597646332), Vector[column](%id = 36893490386597647052), Vector[column](%id = 36893490386597647772), Vector[column](%id = 36893490386597648252), Vector[column](%id = 36893490386597648972), Vector[column](%id = 36893490386597649452), Vector[column](%id = 36893490386597649932), Vector[column](%id = 36893490386597650172)

Although we constructed these cuts relative to a specific spanning tree, they can be shown to be all the cuts of the graph. Consider the set of all edges of the graph. It might seem that this is a cut that has been missed, because it separates the graph into two or more (actually 4) components. However, looking at the edges


{{1, 2}, {1, 3}, {1, 4}, {2, 3}, {3, 4}}

and recalling the definition of a cut, we see that this is not actually a cut. We need the edges in a cut all to have one end in one set of vertices and the other end in the set of other vertices, and this is not possible here. (Maple's IsCutSet is not this strict and returns true.) We can show it is not possible with IsBipartite.



Although we have found all cuts, we now need to test them to find which ones are minimal. The first case with no edges does not cut the graph into two components and we reject it. (It is usually included as a cut to make the cuts into a vector space.)  The 6th one is not minimal. Its first and third entries are 1, meaning it is the ring sum of the first and third fundamental cuts. It splits the graph into three components:

cut := `⊕`(f[1], f[3]); DrawGraph(DeleteEdge(G, cut, inplace = false), size = [200, 200])

{{1, 2}, {1, 4}, {2, 3}, {3, 4}}

This is a case where the ring sum is a disjoint union (not minimal), and arises because f[1] and f[3] have no edges in common:

`intersect`(f[1], f[3])


For the present graph, all but the empty set and `⊕`(f[1], f[3]) are minimal, and so there are 6 minimal cuts. In vector form:

cutvecs := convert(map(tovec, mincuts), set)[]

Vector[column](%id = 36893490386597690532), Vector[column](%id = 36893490386597690652), Vector[column](%id = 36893490386597690772), Vector[column](%id = 36893490386597690892), Vector[column](%id = 36893490386597691012), Vector[column](%id = 36893490386597691132)

So the algorithm is to generate all cuts and test them for minimality. If we neglect the empty cut, then we have to generate 2^r-1 = 2^(n-1)-1 cuts for testing. This is the same as the number of partitions that we generated in the partition algorithm. However, for the fundamental cuts we do not need to test that the components are connected, and more importantly, the ring sum is a more efficient way of finding the edges than partitioning and then using InducedSubgraphs.


The tests for minimality can be made more efficient by finding cases that can be classified without having to check the number of components. Two simple cases with modest savings (implemented here) are:

1. The fundamental cuts are minimal by construction.

2. If two cuts have no edges in common, then the ring sum is a disjoint union and is not minimal. This is the case for example for f[1] and f[3].


There may be other ways to classify without doing the more difficult minimality test, for example variations of (2.) for ring sums with more than two cuts. Hovever, the cutset algorithm is already significantly more efficient than the partition algorithm. Perhaps working with the vectors would be more efficient. Other more efficient algorithms may be known that I am not aware of.


The MinEdgeCuts algorithm in the startup code is reproduced here:


  local mincutsets,n,edges,treeedges,treeG,
  uses GraphTheory;
  if not IsConnected(G) then return Vector([]) end if;
  # choose a spanning tree and find corresponding fundamental cutsets
  r:=n-1;    # =nops(treeedges) = GraphRank(G)
  mincutsets:=table();  # first r ones are the fundamental ones
  for j to r do
    mincutsets[j]:=edges minus partitionedges;
  end do;
  # now find ringsums of all subsets of the set of fundamental cutsets
  fsets := Iterator:-BinaryGrayCode(r,'rank'=2); #skip empty set
  for ref in fsets do;
    if nf=1 then              
      next      # fundamental cutsets already in mincutsets
    elif nf=2 and (fcutsets[1] intersect fcutsets[2] = {}) then
      next      # pair of disjoint cutsets not minimal
    # todo other detections of non-minimal cutsets
    else        # rest the hard way
      for i from 2 to nf do
              cutsetedges:=(cutsetedges union fc) minus (cutsetedges intersect fc);
      end do;
      if nops(vertexpartition)=2 then # cutset is minimal
      end if
    end if;
  end do;
end proc:


Download CutSets.mw



Recently @Nour asked (in part) if Maple can simplify 2*(x^2 - x*y - x*z + y^2) to (x - y)^2 + (y - z)^2 + (z - x)^2. User @vv pointed out that a sum of squares form is not unique. I'm not sure if @Nour's intent was to show that this was always non-negative, but even if the sum of squares form is not unique, it does help to demonstrate that the expression will be nonnegative for all real x,y,z. This reminded me that I have over the years had cases where I wanted to know if a quadratic form (multivariate polynomial with all terms of degree 2) would always be non-negative (a quadratic form will be zero if all variables are zero, so we know it won't always be positive). I have typically done this by hand in the following way, which I demonstrate for the simpler case



We make a (column) vector of the variables, v = [x,y] and then construct the symmetric Matrix A for which Transpose(v).A.v is the required expression. For the squared term 2*x^2, the coefficient 2 is put in the diagonal entry A[1,1] and likewise 4*y^2 leads to A[2,2] = 4. The term -3*x*y comes from the sum A[1,2]+A[2,1] and so half of the coefficient -3 goes to each of A[2,1] and A[1,2]. The Matrix is therefore


A := Matrix(2, 2, {(1, 1) = 2, (1, 2) = -3/2, (2, 1) = -3/2, (2, 2) = 4})


The Matrix is found to be positive definite (positive eigenvalues), which means that the quadratic form is positive for all non-zero real vectors v.



Vector[column]([[4.802775638], [1.197224362]])

One way to get the sum of squares form is to use the Cholesky decomposition (for positive definite or semidefinite matrices), which yields A = L.Transpose(L) = L.I.Transpose(L), and then note that Transpose(L).v is a change of variables to a diagonal quadratic form


L := Matrix(2, 2, {(1, 1) = 2^(1/2), (2, 1) = -(3/4)*2^(1/2), (2, 2) = (1/4)*46^(1/2)}, storage = triangular[lower], shape = [triangular[lower]])

Matrix([[2, -3/2], [-3/2, 4]])

so we have


v2 := Vector(2, {(1) = 2^(1/2)*x-(3/4)*2^(1/2)*y, (2) = (1/4)*46^(1/2)*y})



and we have rewritten the quadratic form in an explicitly nonnegative form.


Another way to the sum of squares would be via the othogonal similarity A = U.Lambda.Transpose(U), where U is an orthogonal matrix of eigenvectors and Lambda is diagonal with eigenvalues on the diagonal. This form works even if the matrix is not positive (semi)definite.


In this simple case, Maple has some other ways of showing that the quadratic form is non-negative for real x,y



Without assumptions, is correctly decides p1 is not always nonnegative, and an example of this for complex x and y is easily found. With assumptions, is fails.

is(p1,nonnegative) assuming real;




At least in Maple 2017, solve suggests any x and y will do, which is not quite right


{x = x, y = y}

In the next worksheet, I implement this method in a general procedure IsQformNonNeg, which works as well as IsDefinite does. For most cases of interest (integer or rational coefficients) there are no problems. For floating point coefficients, as always, numerical roundoff issues arise, and I dealt with this by giving warnings. (Rank is more robust in these cases than IsDefinite, but since this case can be better dealt with by converting floats to rationals, I haven't worked too hard on this - perhaps IsDefinite will be improved.)


The original example 2*(x^2 - x*y - x*z + y^2)  is dealt with as an example in the next worksheet



Download QformPost.mw

@ianmccr posted here about making help for a package using makehelp. Here I show how to do this with the HelpTools package.

The attached worksheet shows how to create the help database for the Orbitals package available at the Applications Center or in the Maple Cloud. The help pages were created as worksheets - start using an existing help page as a model - use View/Open Page As Worksheet and then save from there. The topics and other information are entered by adding Attributes under File/Document Properties - for example for the realY help page these are:

Active=false means a regular help page; Active=true means an example worksheet.

There may be several aliases, for example the cartesion help page also describes the fullcartesian command and so Alias is: Orbitals[fullcartesian],cartesian,fullcartesian and Keyword is: Orbitals,cartesian,fullcartesian

Once a worksheet is created for each help page they are assembled into the help database with the attached file. More information is in the attached file


Page 1 of 1