dharr

Dr. David Harrington

3813 Reputation

17 Badges

18 years, 153 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 replies submitted by dharr

@Thomas Dean 

restart

with(GroupTheory)

Matrix has to have the entries coded as 1..4

v := [1, -1, I, -I]; M := Matrix(4, 4, proc (i, j) options operator, arrow; eval(v[i]*v[j], `~`[`=`](v, [`$`(1 .. 4)])) end proc)

v := [1, -1, I, -I]

Matrix(%id = 36893490351220485524)

G := CayleyTableGroup(M)

_m2203734088608

CayleyTable(G)

Array(%id = 36893490351220474092)

lbl := proc (i) options operator, arrow; v[i] end proc

proc (i) options operator, arrow; v[i] end proc

DrawCayleyTable(G, labels = lbl)

NULL

Download Group.mw

I tried to find this out a while ago, but concluded it wasn't possible. You can use FileTools:-ListDirectory to find a list of all *.mw files in the interface(worksheetdir) directory, but there seems not to be a simple way to find which one is running. I suppose there must be some system-dependent system calls to figure out which programs are running.

@mmcdara Much, much simpler than mine!

@gharouce So, to be clear, you never want the actual numbers out, but just to evaluate some special algebra in which e^2=1, e^3=e etc for one specific variable e in an expression? That is easy to do if your expression is a polynomial in e, but requires more work if you want it to work for, say, e/sin(e^2). Here's the polynomial case:

restart

evalpm:=proc(expr,var)
  local p,sum,power,term,terms;
  p:=collect(expr,var);
  if not type(p,polynom(anything,var)) then error "%1 not polynomial in %2",expr,var end if;
  if nops(p)=1 then terms:=[p] else terms:=[op(p)] end if;
  sum:=0;
  for term in terms do
    power:=diff(term,var)/term*var;
    sum:=sum+ifelse(power::even,coeff(term,var,power),var*coeff(term,var,power));
  end do;
  sum;
end proc:

q := X^3*Y+X^2*Z+X; q2 := (e+1)^3

X^3*Y+X^2*Z+X

(e+1)^3

evalpm(q2, e)

4+4*e

evalpm(q, X)

X*Y+X+Z

NULL

Download Xn.mw

For the more general case

evalpm:=proc(expr,var)
        local p,pwr;
        p:=collect(expr,var);
        evalindets(p,`^`,
            proc(pwr) local res; 
              if op(1,pwr)=var then
                if op(2,pwr)::even then
                    res:=1
                elif op(2,pwr)::odd then
                    res:=var
                else
                    res:=pwr
                end if
              else
                res:=pwr;
              end if;
               res;
             end proc
        );
   end proc:

 

@gharouce You can only set X to one value at a time. But you can do something like this:

s:=X^3*Y + X^2*Z + X;
eval(s,X=1);
eval(s,X=-1);

or you can do something like this:

solve({s=X^3*Y + X^2*Z + X,X^2=1},{s,X});

which gives the two possibilities

@Ronan I originally had ModuleUnload remove the types, thinking that when I did unwith() the module would unload and the types would be removed. But the types still existed - unwith only gets rid of the global names (exports) the package created. The ModuleUnload procedure "is called when the module is destroyed. A module is destroyed either when it is no longer accessible and is garbage collected, or when Maple exits".

So my conclusion is that it isn't really useful.

@Ronan Here's a working example:

restart

packagedir:=interface(worksheetdir);

"C:\Users\dharr\Desktop"

MyThings:=module()
  option package;
  export mysqr,mysqrt;
  local ModuleLoad;
  uses TT=TypeTools;

  ModuleLoad:=proc()
   TT:-AddType(mytype,{identical(x),posint});
  end proc;

  mysqr:=proc(x::mytype);
    x^2
  end proc;

  mysqrt:=proc(x::mytype);
    sqrt(2)
  end proc;
 
end module:
 

libfile:=cat(packagedir, "/MyThings.mla");
march( 'create', libfile, 20 );
savelib('MyThings',libfile);

"C:\Users\dharr\Desktop/MyThings.mla"

 

NULL

Download MyThings.mw

 

restart;

packagedir:=interface(worksheetdir);

"C:\Users\dharr\Desktop"

libname:=libname,packagedir:

with(MyThings);

[mysqr, mysqrt]

mysqr(y);

Error, invalid input: MyThings:-mysqr expects its 1st argument, x, to be of type mytype, but received y

mysqrt(2.5);

Error, invalid input: MyThings:-mysqrt expects its 1st argument, x, to be of type mytype, but received 2.5

mysqr(x);

x^2

NULL

 

Download test.mw

Deciding whether or not a given cut is minimal (breaks the graph only into two components) is the slow part for the above algorithm; the test uses ConnectedComponents that needs to perform an algorithm that searches through the graph. In the description of the above algorithm two cases were mentioned where this hard test could be avoided. One was that the fundamental cuts were minimal by construction. The second was for ring sums between two fundamental cuts: if the ring sum was between two edge sets with no common edges, it was not minimal and the hard test was not required.

I speculated that this might be possible to simplify for ringsums of more than two cuts. The following algorithm accumulates the ringsums incrementally, and at each stage tests for lack of common edges.

We need to test all 2^r subsets of the set of r fundamental cuts (except the empty set). Let's code the subsets by a list of integers that index the fundamental sets. For example the index [1,2,4,5] will mean evaluating the ring sum of the first, second, fourth and fifth subsets. The subsets iterator in the combinat package can iterate through all subsets of integers in the list [1,2,3,...,r], thereby successively finding all subsets. More importantly, it does it in an order that means we can use previously evaluated ringsums to find new ones. The order is all subsets with zero elements, then one element, then two elements, etc., and in numerical order within those cases, e.g., for r=3 this is [], [1], [2], [3], [1,2], [1,3], [2,3], [1,2,3]. So as an example (for r>4), when we come to evaluate the ringsum of subsets [1,2,3,7,9] we find the ringsum of subsets [1,2,3,7] has already been done, and so we can evaluate the ringsum between the [1,2,3,7] result and fundamental cut [9], and store the result in cutsets[1,2,3,7,9] for later use. If the [1,2,3,7] and [9] edge sets have no common edges, this will not be minimal and we can avoid the hard test in this case.

This was implemented using a Maple table cutsets that was indexed with the lists of integers. All 2^r cutsets need to be stored, so there is a larger memory requirement. Profiling the algorithm on @Icz's graph with 20 vertices  and 72 edges and 314415 minimal cuts showed that the time saved by incrementally doing the ringsums was approximately offset by the extra complexity of the algorithm. For this graph at least, although many more hard tests were avoided, this was still a small fraction of the total.

                       Old         New
fundamental             19          19
easy (no common edges) 104       25063
hard tests          524164      499205
total = 2^19-1      524287      524287

The total times were appromimately the same. 342 out of the 370 seconds (92%) was spent on hard tests, so those are still limiting here. The number of hard tests was reduced by only 5%.

Edit: This analysis suggests that testing all pairs out of 2^r may be the way to go, or perhaps avoid a brute-force method using the much more complicated algorithm of @Icz's ref [1].

Update: testing all pairs is much, much slower, so the mixed strategy here or above seems to be optimum for a brute force method. 

restart

MinimalEdgeCuts:=proc(G1::GRAPHLN)
  local mincutsets,edges,treeedges,G,treeG,
    i,j,r,vertexpartition,partitionedges,nf,fc,
    fsets,ref,fcutsets,cutsetedges,cutsets,iter,f1,f2,f1intf2,f1ringf2,s;
  uses GraphTheory;
  if not IsConnected(G1) then error "graph is not connected" end if;
  G:=UnderlyingGraph(G1); # make undirected, unweighted, without loops
  r:=NumberOfVertices(G)-1;
  edges:=Edges(G);
  cutsets:=table();     # all 2^r cuts
  mincutsets:=table();  # first r ones are the fundamental ones
  iter:=combinat:-subsets([$1..r]); #iterate through all 2^r
  iter[nextvalue](); #skip empty list
  # choose a spanning tree
  treeG:=SpanningTree(G);
  treeedges:=Edges(treeG);
  # find r corresponding fundamental cutsets, one for each tree edge
  for j to r do
    iter[nextvalue]();  
    vertexpartition:=ConnectedComponents(DeleteEdge(treeG,treeedges[j],'inplace'=false));
    partitionedges:=`union`(map(Edges,map2(InducedSubgraph,G,vertexpartition))[]);
    mincutsets[j]:=edges minus partitionedges;
    cutsets[j]:=mincutsets[j];
  end do;
  j:=r;
  # now find ringsums of all subsets of the set of fundamental cutsets
  while not iter[finished] do
    s:=iter[nextvalue]();
    f1:=cutsets[op(1..-2,s)];
    f2:=cutsets[s[-1]];
    f1intf2:=f1 intersect f2;
    f1ringf2:=(f1 union f2) minus f1intf2;
    cutsets[op(s)]:=f1ringf2;
    if f1intf2 = {} then
      next
    else
      # hard test
      vertexpartition:=ConnectedComponents(DeleteEdge(G,f1ringf2,inplace=false));
      if nops(vertexpartition)=2 then # cutset is minimal
        j:=j+1;
        mincutsets[j]:=f1ringf2
      end if
    end if;
  end do;
  Vector(j,mincutsets);
end proc:

with(CodeTools:-Profiling):

Profile(MinimalEdgeCuts);

g := "S~tIID@OI?{@n~V?goYEDOWd?qI?sJ?[C";
G := GraphTheory:-ConvertGraph(g);
#GraphTheory:-DrawGraph(G);

"S~tIID@OI?{@n~V?goYEDOWd?qI?sJ?[C"

GRAPHLN(undirected, unweighted, [1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13, 14, 15, 16, 17, 18, 19, 20], Array(1..20, {(1) = {2, 3, 4, 5, 12, 13}, (2) = {1, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 14, 15, 16, 17, 18, 19, 20}, (3) = {1, 2, 4, 12, 13, 20}, (4) = {1, 2, 3, 5, 13, 14}, (5) = {1, 2, 4, 6, 13, 14}, (6) = {2, 5, 7, 13, 14, 15}, (7) = {2, 6, 8, 13, 15, 16}, (8) = {2, 7, 9, 13, 16, 17}, (9) = {2, 8, 10, 13, 17, 18}, (10) = {2, 9, 11, 13, 18, 19}, (11) = {2, 10, 12, 13, 19, 20}, (12) = {1, 2, 3, 11, 13, 20}, (13) = {1, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 14, 15, 16, 17, 18, 19, 20}, (14) = {2, 4, 5, 6, 13, 15}, (15) = {2, 6, 7, 13, 14, 16}, (16) = {2, 7, 8, 13, 15, 17}, (17) = {2, 8, 9, 13, 16, 18}, (18) = {2, 9, 10, 13, 17, 19}, (19) = {2, 10, 11, 13, 18, 20}, (20) = {2, 3, 11, 12, 13, 19}}), `GRAPHLN/table/1`, 0)

mincutsets:=CodeTools:-Usage(MinimalEdgeCuts(G)):
numelems(mincutsets);

memory used=22.76GiB, alloc change=447.40MiB, cpu time=6.16m, real time=3.88m, gc time=2.83m

314415

PrintProfiles(MinimalEdgeCuts);

MinimalEdgeCuts
MinimalEdgeCuts := proc(G1::GRAPHLN)
local mincutsets, edges, treeedges, G, treeG, i, j, r, vertexpartition,
  partitionedges, nf, fc, fsets, ref, fcutsets, cutsetedges, cutsets, iter, f1
  , f2, f1intf2, f1ringf2, s;
     |Calls Seconds  Words|
PROC |    1 369.781 3054875702|
   1 |    1   0.000    520| if not GraphTheory:-IsConnected(G1) then
   2 |    0   0.000      0|     error "graph is not connected"
                            end if;
   3 |    1   0.000  90679| G := GraphTheory:-UnderlyingGraph(G1);
   4 |    1   0.000     10| r := GraphTheory:-NumberOfVertices(G)-1;
   5 |    1   0.000    837| edges := GraphTheory:-Edges(G);
   6 |    1   0.000    263| cutsets := table();
   7 |    1   0.000    263| mincutsets := table();
   8 |    1   0.000   3528| iter := combinat:-subsets([$ 1 .. r]);
   9 |    1   0.000     72| iter[nextvalue]();
  10 |    1   0.016   2570| treeG := GraphTheory:-SpanningTree(G);
  11 |    1   0.000    192| treeedges := GraphTheory:-Edges(treeG);
  12 |    1   0.000      0| for j to r do
  13 |   19   0.000   1432|     iter[nextvalue]();
  14 |   19   0.000  22390|     vertexpartition := GraphTheory:-
                                  ConnectedComponents(GraphTheory:-DeleteEdge(
                                  treeG,treeedges[j],('inplace') = false));
  15 |   19   0.000  24057|     partitionedges := `union`(map(GraphTheory:-
                                  Edges,map2(GraphTheory:-InducedSubgraph,G,
                                  vertexpartition))[]);
  16 |   19   0.000   1710|     mincutsets[j] := edges minus partitionedges;
  17 |   19   0.000    266|     cutsets[j] := mincutsets[j]
                            end do;
  18 |    1   0.000      0| j := r;
  19 |    1   3.621      4| while not iter[finished] do
  20 |524268  10.529 58981573|     s := iter[nextvalue]();
  21 |524268   1.023 6028911|     f1 := cutsets[op(1 .. -2,s)];
  22 |524268   0.548 1048536|     f2 := cutsets[s[-1]];
  23 |524268   3.079 9436865|     f1intf2 := f1 intersect f2;
  24 |524268   4.557 50878477|     f1ringf2 := f1 union f2 minus f1intf2;
  25 |524268   0.875 4414435|     cutsets[op(s)] := f1ringf2;
  26 |524268   1.793 1572856|     if f1intf2 = {} then
  27 |25063   0.000      0|         next
                                else
  28 |499205 341.972 2917786560|         vertexpartition := GraphTheory:-
                                      ConnectedComponents(GraphTheory:-
                                      DeleteEdge(G,f1ringf2,inplace = false));
  29 |499205   1.013 1497615|         if nops(vertexpartition) = 2 then
  30 |314396   0.267      0|             j := j+1;
  31 |314396   0.473 2762501|             mincutsets[j] := f1ringf2
                                    end if
                                end if
                            end do;
  32 |    1   0.015 318580| Vector(j,mincutsets)
end proc
 

NULL

Download MinEdgeCutsNewIterProfile.mw

@lcz The new method discussed here is about 35% faster for your graph.

@lcz I'm not sure how to interpret the claims in Ref [1]. It doesn't have many citations, and a quick look at those suggests they are not about improved algorithms, so perhaps it is not a significant paper (or perhaps it is just unnoticed).

I improved the efficiency slightly in a new version. There are 2^(n-1) partitions to test for minimality, and IsConnected() does a breadth-first search, so will depend on size also. 

I am looking at doing this through combining fundamental cutsets. There is an algorithm to find these from a spanning tree, and combine them to make all (not just minimal) 2^(n-1) cutsets. Then there will be the same number of minimality tests, but it might be possible to easily do some, or the test may be more efficient.

@lcz The cutset problem between two vertices seems more common, since it has applications in fault theory in networks, which made it hard to find references for the above problem.

@dharr The Print routine for Iterator:-Combination incorrectly shows zero-relative ranks, which led to a duplcate partition in my original code. (SCR submitted.)

restart

M := Iterator:-Combination(4, 2)

_m1640374251104

Following suggests rank starts at 0

Print(M, showrank)

0: 0 1
1: 0 2
2: 1 2
3: 0 3
4: 1 3
5: 2 3

But actually it starts at 1

for i in M do
  Rank(M),i;
end do;

1, Vector[row](2, {(1) = 0, (2) = 1})

2, Vector[row](2, {(1) = 0, (2) = 2})

3, Vector[row](2, {(1) = 1, (2) = 2})

4, Vector[row](2, {(1) = 0, (2) = 3})

5, Vector[row](2, {(1) = 1, (2) = 3})

6, Array(%id = 36893489787840041436)

So to get only the second half, start with rank 4, not 3.

NULL

Download IteratorRank.mw

@Oliveira Here it is with the syntax error fixed.

TravelingSalesperson.mw

@PaulNewton You launcher worksheet idea is possible. See ?DocumentTools,RunWorksheet. Then the code following a successful "return" from the launched worksheet does not execute. The launched worksheet runs "headless", meaning that some of its interface functions are not running, so it may not do what you want. 

@mmcdara The original question seemed more about 1-D vs 2-D, so I wasn't thinking much about efficiency. There has been some discussion on how to speed up cartesian products, so I'm not surprised that it is slower:

https://www.mapleprimes.com/posts/41838-Cartesian-Products-Of-Lists

1 2 3 4 5 6 7 Last Page 1 of 36