Carl Love

## 27187 Reputation

11 years, 340 days
Himself
Wayland, Massachusetts, United States
My name was formerly Carl Devore.

## Use the order as an index...

You just need to use the given ordering as an index of the vertex list:

R:= IsChordal(s, eliminationordering):
if [R][1] then op(3,s)[R[2]] fi;

## GAMMA...

If by Gamma you mean the generalized factorial function, then you need to spell it GAMMA (all uppercase). The cause of the error is that it doesn't know what Gamma means,

## Better proc; Symmetry analysis...

I improved my algorithm from above (made it faster) and added extensive explanation of how it works. Considering the code you've posted recently, you desperately need to study this explanation.

Your example graph has a lot of symmetry, 48 automorphisms. In the worksheet below, I show some techniques for analyzing the symmetry, including dividing the 300 minimal dominating sets into 12 equivalence classes under the automorphisms. The equivalence relation is specifically that two minimal dominating sets are equivalent if there is an automorphism of mapping one to the other.

Minimal dominating sets of a graph

 > restart :
 > #Author: Carl Love , the Glorious 1st of June, 2023 # #Returns all minimal dominating sets of the input simple graph. #                                          #Reference #s for expository comments: MinDominatingSets:= proc(G::Graph, {any::truefalse:= false})                # 1 local                              n:= nops(op(3,G)), V:= {\$n},                                            # 2     A:= [seq](op(4,G) union~ `{}`~(V)),                                     # 3     N:= table([()= {}]), N0,                                                # 4     R:= <()>,                                                               # 5     s, Ns, v                                                                # 6 ;     do                                                                      # 7         N0:= eval(N); N:= evaln(N);                                         # 8         for s,Ns in eval(N0) do                                             # 9             for v from `if`(s::posint, s, `if`(s=(), 0, s[-1])) + 1 to n do #10                 if (N[s,v]:= Ns union A[v]) = V then                        #11                     if any then return {s,v} else R,= {s,v} fi              #12                 fi                                 od         od;         if numelems(R) <> 0 then return {seq}(R) fi                         #13     od     end proc :

Expository comments on the code and algorithm:

#0: Note that the procedure does not call any command from GraphTheory. There's usually no need to do that. Indeed, this procedure doesn't use any commands from any package.

#1: If the keyword option any is passed, then the procedure returns only the first minimal dominating set that it finds.

#2: One key to efficient graph theory algorithms is to never refer to vertices by their labels. Instead, they are always referred to by their index into the list of vertices. The list of vertices is op(3, G), the same as GraphTheory:-Vertices(G). So, V here is simply the set {\$1..n} where n is the number of vertices. All further references to vertices, including within edges, will use these indices.

#3: The second key to efficient graph theory algorithms (this is very important) is what I call the adjacency array. I don't know if it has a formal name. It's returned by op(4, G). It's a sparse representation of the adjacency matrix. It's a one-dimensional array A indexed by the vertices such that A[k] is the set of vertices adjacent to k.

For the purpose of calculating dominating sets, it's convenient to consider vertices as being adjacent to themselves; so, the union~ adds each k to the set A[k].

#4: N is a table that is the key data structure of this algorithm. The kth iteration of the main loop (#7) constructs all k-susbsets of V and their neighborhoods. The k-subsets (stored as sequences) are the indices of  N, and the neighborhoods are the entries. They are constructed from the (k-1)-subsets and their neighborhoods, which are moved to table N0 at at the beginning of loop (#7).

#5: R is a array to hold the dominating sets to be returned.

#6: On the kth iteration of the main loop (#7), s will be a (k-1)-subset of vertices as a sequence, Ns will the neighborhood of s, and v be a vertex to be added to s to make a k-subset.

#7: Everything that needs to be said about this line is covered in #4 and #6 above. Although I've referred to the kth iteration of this loop several times, the code itself does not use k.

#8: Move the (k-1)-subsets and their neighborhoods to table N0, then clear and garbage-collect N in preparation for the k-subsets.

#9: When a for loop has two control variables, the first is the index into the container, and the second is the container's element at that index. So, in this case, s runs over all indices of N0, and Ns = N0[s].

#10: Add to s all vertices v in such a way that each subset of V is constructed at most once. Specifically, we only add vertices v that are greater than all vertices is s.

#11: The index of N is the new k-subset (stored as a sequence); the entry is its neighborhood. If the neighborhood is all of V, then the index is a dominating set.

#12: If we have a dominating set and keyword any is set, then return it; otherwise save the dominating set for later return. The R,= is an efficient way to append a new element to an Array without indexing.

#13: If any dominating sets have been stored in R, then return them. They're necessarily all dominating sets of size k and all minimal.

Note that a return statement can be used anywhere in a procedure. You do not need to exit loops before using return.

 > #Example: A 7-regular graph with 24 vertices # g:= (GT:= GraphTheory):-Graph(     [\$24],     {         {1, 2}, {1, 3}, {1, 4}, {1, 5}, {1, 7}, {1, 9}, {1, 11},         {2, 3}, {2, 5}, {2, 6}, {2, 7}, {2, 12}, {2, 14},         {3, 4}, {3, 7}, {3, 8}, {3, 9}, {3, 16},         {4, 5}, {4, 9}, {4, 10}, {4, 11}, {4, 17},         {5, 6}, {5, 11}, {5, 12}, {5, 19},         {6, 7}, {6, 12}, {6, 13}, {6, 14}, {6, 21},         {7, 8}, {7, 14}, {7, 15},           {8, 9}, {8, 14}, {8, 15}, {8, 16}, {8, 22},         {9, 10}, {9, 16}, {9, 17},         {10, 11}, {10, 17}, {10, 18}, {10, 19}, {10, 23},         {11, 12}, {11, 18}, {11, 19},  {12, 13}, {12, 19}, {12, 20},         {13, 14}, {13, 19}, {13, 20}, {13, 21}, {13, 24},         {14, 15}, {14, 21},  {15, 16}, {15, 21}, {15, 22}, {15, 24},         {16, 17}, {16, 22}, {16, 23},  {17, 18}, {17, 22}, {17, 23},         {18, 19}, {18, 20}, {18, 23}, {18, 24},  {19, 20},         {20, 21}, {20, 23}, {20, 24},  {21, 22}, {21, 24},         {22, 23}, {22, 24}, {23, 24}     } );

 > #Get 1 minimal dominating set: ds:= CodeTools:-Usage(MinDominatingSets(g, any));

memory used=1.01MiB, alloc change=0 bytes, cpu time=0ns, real time=13.00ms, gc time=0ns

 > #Get them all: ds:= CodeTools:-Usage(MinDominatingSets(g));

memory used=13.89MiB, alloc change=37.00MiB, cpu time=110.00ms, real time=80.00ms, gc time=62.50ms

 > nops(%);

 > #Check symmetries (automorphisms): # GpT:= GroupTheory: P:= CodeTools:-Usage(     curry~(GpT:-PermApply, [GpT:-Elements(GT:-AutomorphismGroup(g))[]]) ): nops(P);

memory used=2.21MiB, alloc change=0 bytes, cpu time=187.00ms, real time=440.00ms, gc time=0ns

 > #There are 48 automorphisms; the graph has many symmetries. #Partition the 300 dominating sets into #equivalence classes under automorphism:  # EC:= CodeTools:-Usage(     [do ec:= {(`{}`~@op@P~)(ds[1])[]} until (ds minus= ec) = {}] );

memory used=1.12MiB, alloc change=0 bytes, cpu time=31.00ms, real time=21.00ms, gc time=0ns

 > #Check global symmetries: #    1. Given any pair [u,v] of vertices, is there an automorphism of g #    that takes u to v? # evalb(nops((`{}`@op@P)~({\$24})) = 1);

 > #    Yes, there is. #
 > #    2. Given any pair [u,v] of vertices, is there an automorphism of g #    that exchanges them? # pairs:= combinat:-choose([\$24], 2): P2:= subsindets((P~)~(pairs), listlist(posint), `[]`~@op): andseq(pairs[k][[2,1]] in P2[k], k= 1..nops(pairs));

 > #    There is not. #
 > #    3. Given any 2 pairs {u1,v1}, {u2,v2} of vertices, is there an #    automorphism of g taking the 1st pair to the 2nd? # pairs:= combinat:-choose({\$24}, 2): evalb(nops(     subsindets(subsindets((P~)~(pairs), set(list), `{}`~@op), list, `{}`@op) ) = 1);

 > #    There is not. #
 > #    4. Given any two edges of g, is there an automorphism taking the 1st #    to the 2nd? # edges:= GT:-Edges(g): evalb(nops(     subsindets(subsindets((P~)~(edges), set(list), `{}`~@op), list, `{}`@op) ) = 1);

 > #    There is not. #
 > #Classify the edges by their number of automorphisms: # ListTools:-Classify(     nops@`{}`@op@`{}`~@op@P~, edges );

 > nops~(%);

 > #So there are 12 edges with 12 automorphic images and 72 with 24 images.
 >

## Substitute reciprocals...

It can be done simply with subs by substituting the reciprocal of the desired denominator for the reciprocal of the current denominator, like this:

subs(1/rhs(eq_Dx)= 1/Dx, eq_K1_m4);

## Array...

Let's say that you want the indices from 1 to 2, 1 to 3, 1 to 4, and 1 to 5. Then do

S:= Array(1..2, 1..3, 1..4, 1..5);

## Maple Calculator...

I'm guessing that @mmcdara is joking with his "non-answer", but his Reply got me thinking, and I came up with another "non-answer" that may be useful for you (so, it's actually an Answer).

If you have a device (such as a phone) with a camera and the Android or iOS operating system, then download the free Maple Calculator app to it. The app can translate photos of typeset or handwritten mathematical formulae into Maple code. I'm not aware of that photo-translation capability being available in Maple itself.

Tap the camera 📷 icon, then frame and photograph your formula. After a few seconds, you'll probably get an error message like "Sorry, this took too long to compute", which doesn't matter because you're not interested in computing it anyway. Tap the "upload to Cloud" icon. When it's done, it'll tell you the filename in the MapleCloud. In this case, the filename is the exact Maple sum command that you need.

Of course, the intended purpose of this photo-to-cloud facility is to let you download that Cloud file to Maple proper. But, for your immediate purpose, the file title has everything you need, so you don't even need to open Maple.

## Needs to be in a procedure or loop...

The error is because the target label is not in a procedure or loop. Not being in a procedure or loop is also called being "at top level".

There is no good reason for you to ever use goto. In addition to being widely frowned upon in most languages, the Maple implementation is very inefficient.

There used to be a help page for it. But using the command is so discouraged that the help page has been removed.

## evala only helps with radicals and RootO...

As far as I know, evala (and all of its large number of subcommands) only makes any meaningful change to exact expressions containing radicals, explicit rational fractional exponents, and RootOfs of polynomials (not any other RootOf).

## Just with simplify...

The job can be done---obtaining the same result as Tom---without subs and without needing to rearrange your substitution expression---by using a single command: a form of simplify known as "with side relations" (see help page ?simplify,siderels). The exact command in this case is

Note that the 2nd argument must be given as a set, and it can include multiple equations. The 3rd argument is often not needed, but in this case it's needed if you want to obtain your desired expression rather than some other substitution/simplification. If given, it should be a list of variables to eliminate, given in the priority order that you want them eliminated.

## kernelopts(numcpus= 1)...

Set

kernelopts(numcpus= 1)

Many years ago, I discovered a bug such that one could only change kernelopts(numcpus) immediately after a restart. If it was changed afterwards (I guess after the first garbage collection occurs), the kernelopts command would claim that it was changed, but in reality it wasn't. I don't know if this has been fixed.

You wrote:

• My code does not do any mutlithreading. So I do not need it.

You don't know that! You do very high-level stuff that calls thousands of library procedures. Some of those probably attempt multithreading if they see CPUs available for it. But, I've never seen any Maple library command error out due to the absence of multithreading.

## A procedure for it...

Here's a procedure for it. It returns a domination set of minimal size so that its dominance can be easily verified. Then it's trivial to count it if you want.

 > restart :
 > #Author: Carl Love , the Glorious 1st of June, 2023 # DominationSet:= (G::Graph)-> local V:= {\$nops(op(3,G))}, A:= op(4,G) union~ `{}`~(V), s, v, U:= table([{}= {}]), Us;     in V do         for s in {indices}(U, 'nolist') do             Us:= U[s]; U[s]:= evaln(U[s]);             for v in V minus s do                 if (U[s union {v}]:= Us union A[v]) = V then return {s[], v} fi             od         od     od :
 > #Example: ed:={{1, 2}, {1, 3}, {1, 4}, {1, 5}, {1, 7}, {1, 9}, {1, 11}, {2, 3}, {2, 5}, {2, 6}, {2, 7}, {2, 12}, {2, 14},  {3, 4}, {3, 7}, {3, 8}, {3, 9}, {3, 16}, {4, 5}, {4, 9}, {4, 10}, {4, 11}, {4, 17}, {5, 6}, {5, 11}, {5, 12}, {5, 19}, {6, 7}, {6, 12}, {6, 13}, {6, 14}, {6, 21}, {7, 8}, {7, 14}, {7, 15}, {8, 9}, {8, 14}, {8, 15}, {8, 16}, {8, 22}, {9, 10},  {9, 16}, {9, 17}, {10, 11}, {10, 17}, {10, 18}, {10, 19}, {10, 23}, {11, 12}, {11, 18}, {11, 19}, {12, 13}, {12, 19}, {12, 20}, {13, 14}, {13, 19}, {13, 20}, {13, 21}, {13, 24}, {14, 15}, {14, 21}, {15, 16}, {15, 21}, {15, 22}, {15, 24}, {16, 17}, {16, 22}, {16, 23}, {17, 18}, {17, 22}, {17, 23}, {18, 19}, {18, 20}, {18, 23}, {18, 24}, {19, 20}, {20, 21}, {20, 23}, {20, 24}, {21, 22}, {21, 24}, {22, 23}, {22, 24}, {23, 24}}: g:= GraphTheory:-Graph(ed);

 > ds:= CodeTools:-Usage(DominationSet(g));

memory used=3.66MiB, alloc change=0 bytes, cpu time=32.00ms, real time=25.00ms, gc time=0ns

 > op(4,g)[[ds[]]];

 > seq[reduce= `union`](%) union ds; nops(%);

 >

## Rational functions...

You wrote:

• Since I cannot feed non-algebraic expressions into polynomial solvers, I want to extract only their numerators (which should be algebraic)....

You are misusing the mathematical term algebraic. Since I think that there's a good chance that such misuse will lead you to make even-more-serious errors, I'm going to digress for a bit to attempt to correct that misuse. First, note that Maple has a type named algebraic. That type has nothing to do with the mathematical concept! This discussion is strictly about the mathematical concept.

I'll assume that you understand the definitions of all these mathematical objects; let me know if you need clarification on any:

• integer
• rational number
• polynomial
• polynomial in the variables {x, y, etc.}
• polynomial in the variables {x, y, etc.} with integer coefficients
• polynomial in the variables {x, y, etc.} with coefficients that are polynomials in the variables {a, b, etc.} with integer coefficients.

Your terminology usage suggests that the above is perhaps the limit of your knowledge of the precise definitions of even-larger classes of objects. Thus:

Definition: A rational function is a quotient of polynomials, or an expression that can be re-expressed as such via reversible common arithmetic operations (as learned in high-school algebra). As above, we may refine the class by specifying the variables and the type of coefficients.

Definition:  An algebraic function is a rational function with the possible inclusion of non-integer rational-number (fractional) exponents (such as square roots). These exponents may appear not just on individual variables but also on any subexpressions.

Let me know if you need any clarification of those definitions.

Now, returning to your immediate questions, it seems that by "non-algebraic" you actually mean non-polynomial.

• These solutions should also solve the non-algebraic [non-polynomial] system then.

Almost true. The solutions thus obtained may be a nontrivial superset of the solutions of the original system. They need to be verified, which is a much much easier process than obtaining that superset. In the case at hand, your original system is rational functions, and thus the only thing that really needs to be verified is that the solutions do not make any of the original denominators zero.

• To extract the numerators, I am using ((numer@evala@:-Norm@numer)~@eval)(Eqs)....

Obviously, that's either a command that I composed, or a variation of such. In this case, all that you need is numer~(Eqs) because they are rational functions. And the eval is only needed if you want to instantiate parameters. But deconstructing anyway:

1. @ is the function-composition operator: (f@g)(x,y) = f(g(x,y)). As you may recall from a long-ago math course, function composition is associative but commutative. Let me know if you need a refresher.
2. is the elementwise operator (see help page ?elementwise). It causes the operator to which it is appended to be applied individually to the elements (members) of a container (such as a set) that follows. In this case specifically, it applies the operators to your individual equations. (By equation of course I mean an expression implicitly equated to 0.)
3. evala is a "package" (in an older Maple sense of "package") of commands for performing sophisticated operations on algebraic functions. Since you only have rational functions in this case, it's not needed. When I originally gave you the command, the situation was much more complicated with numerous radicals in the equations. Those were algebraic, but not rational, functions.
4. Norm is a command in evala for converting algebraic functions to rational functions. This process necessarily involves irreversible operations (such as squaring both sides of an equation to remove square roots). Of course, it's not needed in this case.
5. :- is a purely syntactic "disambiguation" operator. On the off-chance that you had loaded some other package (such as LinearAlgebra) that also has a command named Norm (totally unrelated), the :- specifies that that is not the Norm that we mean.

• 1. Am I wrong?

No. Other than the terminology issue and using a needlessly complicated command for the easy job of extracting the numerators of rational functions, you're totally correct.

• 3. Why length(((numer@evala@:-Norm@numer)~@eval) numer~(Eqs)) is much larger than length(Eqs)? Shouldn't it be smaller since I just extract the numerators?

Example: Let e:= 1/x + 1/y + 1/z. The length of e is 25. The numerator of e is y*z + x*z + x*y. It's length is 37. Although it's not important at this point to understand how the exact values 25 and 37 are obtained (it's idiosyncratic to Maple), isn't it obvious that the numerator alone should be considered "lengthier" than the original? If it isn't obvious, extend e to more terms and see what happens.

Finally, your title also mentions "simplify". Other than instantiating parameters to numeric values, I wouldn't simplify, unless absolutely necessary due to memory usage. Compared to the solving process, it's trivial for the polynomial solver to simplify in whatever way works best for it.

If you don't instantiate at least some parameters, I think there's very little hope of obtaining any solutions. I'm not saying that those solutions don't exist, just that I think that they're not possible to obtain symbolically.

## F makes no assignment to f...

Your procedure F does not make any assignment. In particular, it makes no assignment to f. Thus, no matter how many times you call F(i), and no matter what value you use for i, the only assigned values of f will be f[2] and f[3]. To correct this, the statement in the for loop should be

f[i]:= F(i)

## Not intended for strings...

If you read the extremely simple and extremely short code of MmaTranslator:-Mma:-LeafCount, then I think that you'll immediately see what's going wrong. Basically, it attempts to break everything down into numbers and names. Your input is a string, so it cannot be broken down into numbers and names. If you remove the quotes from your input, then you'll get the that you expect. If you still can't understand it after reading the code, let me know, and I'll explain it in more detail.

The entire Mma subpackage of MmaTranslator reads to me as if it was written by a rank amateur, such as someone who had 1 or 2 programming courses freshman year at university, then maybe self-studied Maple for a month.

## One plot per file...

There can only be one "plot" per file. I put "plot" in quotes because any number of plots can be combined into a single plot by passing an Array of plots (rather than a sequence of plots) to plots:-display. This can be done simply by enclosing the seq command in angle brackets: <seq(...)>.

Example:

plotsetup('png', 'plotoutput'= "ArrayOfPlots", 'plotoptions'= "width=1920,height=1080"):
cols1:= ["red", "green", "blue"]:
plots:-display(<seq(plot([\$20], 'color'= cols1[j]), j= 1..3)>);
plotsetup('default');
#Reset output to normal screen

I don't know whether the above will be useful to you. If you need to afterwards separate the individual plots, then it may not be very useful. In that case, you can create essentially an "array of files", 1 plot per file, like this:

restart:
cols1:= ["red", "green", "blue"]:
Plots:= <seq(plot([\$20], 'color'= cols1[j]), j= 1..3)>:
PlotToFile:= proc(P::specfunc(PLOT), j::posint)
plotsetup(
'png',
'plotoutput'= cat("ArrayOfPlots_", j),
'plotoptions'= "width=1920,height=1080"
);