dharr

Zeon algebra...

Answered: dharr 8445

June 21 2023

1 0

See also the thread here
How-To-Obtain-All-Cycles--In-A-Graph

>

Using the adjacency matrix and zeons to count cycles.

See R. Schott and G.S. Staples, Complexity of counting cycles using zeons, Comp. Math. Appl. 62 (2011) 1828–1837, doi: 10.1016/j.camwa.2011.06.026

Zeons implemented with bitwise operations on Maple integers in the Bits package.

>

Setting a bit in an integer.

>

Warning, (in SetBits) `i` is implicitly declared local

>

Going back

>

The algorithm starts with two versions of the adjacency matrix. Matrix B has the (i,j) entry as a coded version of i and j that is the beginning of a path or cycle.

Matrix A entry (i,j) has only the coded versionsof j; this is the vertex that will be added to the path at each step by evaluating B.A with a special rule for multiplication.

At any stage the list in entry (i,j) of B contains integers representing possible paths between i and j of a certain length k.

>

Adj := AdjacencyMatrix(G); A := Matrix(n, n, proc (i, j) options operator, arrow; if Adj[i, j] = 1 then SetBit(j) else 0 end if end proc); B := Matrix(n, n, proc (i, j) options operator, arrow; if Adj[i, j] = 1 then [Or(SetBit(i), SetBit(j))] else [] end if end proc); B, A

Matrix(%id = 36893490887404493628), Matrix(%id = 36893490887404491348)

Consider a list of two paths of length 2 from vertex 1 to vertex 3: 1-2-3 and 1-7-3 that would be in entry (1,3). Part of the next step would be trying to add vertex 7 to each of these paths.

>

The multiplication of two entries in the process of Matrix multiplication adds the a vertex to each of the possibilities in b. If it is a repeated vertex, then the possibility disappears. A duplicate is detected as the And of two values being nonzero. Adding the vertex is done by Or

>

And(7, 64); And(69, 64); Or(7, 64); Decode(%)

Note that although we add the vertex, we lose information about the order. So this algorithm cannot output the cycles themselves. The following implements this multiplication.

>	map((x,y)-> if Bits:-And(x,y)=0 then Bits:-Or(x,y) else NULL end if,b,a);

A repeatng vertex isn't necessarily a cycle; only if the starting and ending vertex are the same. This will happen if we are multplying B[i,j]*A[j,p] with i=p, i.e., we are calculating a diagonal entry.

To prevent double counting, we only keep the cycle if the smallest vertex is i. The full multiplication rule is

>

mult:=(x,y)->if y = 0 or x=[] then []
             else
               map((u,w)-> if Bits:-And(u,w)=0 then
                              return Bits:-Or(u,w)
                           else
                              if i=p and i=Bits:-FirstNonzeroBit(u) + 1 then record_cycle end if;
                              return NULL;
                           end if,
                    x,y)
             end if:

Implement the Matrix multiplication with the above routine

>	MM:=proc(B,A) local C,i,j,p; C:=Matrix(n,n): for i to n do # rows of B diving down for p to n do # cols of A C[i,p]:=[seq(mult(B[i,j],A[j,p])[],j=1..n)]; #Row(B,i).Column(A,p) end do end do: C end proc:

>

Matrix(%id = 36893490887341852180)

Next one

>

B := C; C := MM(B, A); Show(C)

Matrix(%id = 36893490887341836396)

And again

>

B := C; C := MM(B, A); Show(C)

Matrix(%id = 36893490887253294844)

So here is a procedure to count cycles in a graph or digraph.

Output is a sparse table, indexed by the number of vertices in the cycles

Use add(eval(CycleCounts(G))) to find the total number of cycles.
This includes 2-cycles, one for each edge, so subtract the number of edges if you do not want 2-cycles in the total.

For undirected graphs, some authors count each cycle twice, one in each orientation, but here each cycle is counted only once.

>

CycleCounts:=proc(G::Graph)
  uses GraphTheory;
  local k,A,B,C,p,n,i,j,mult,cycles,SetBit;
  cycles:=table('sparse');
  n:=NumberOfVertices(G);
  SetBit := m-> Bits:-Join(convert(Vector(n, {m = 1}, 'fill' = 0), list));
  mult:=(x,y)->if y = 0 or x = [] then []
             else
               map((u,w)-> if Bits:-And(u,w)=0 then
                              return Bits:-Or(u,w)
                           else
                              if i=p and i=Bits:-FirstNonzeroBit(u) + 1 then cycles[k]++ end if;
                              return NULL;
                           end if,
                    x,y)
             end if;
  A:=Matrix(n,n, (i, j)-> if AdjacencyMatrix(G)[i, j] = 1 then SetBit(j) else 0 end if);
  B:= Matrix(n,n, (i, j)-> if AdjacencyMatrix(G)[i, j] = 1 then {Bits:-Or(SetBit(i),SetBit(j))} else [] end if);
  for k from 2 to n do
    C:=Matrix(n,n):
    for i to n do # rows of B diving down
      for p to n do # cols of A
        C[i,p]:=[seq(mult(B[i,j],A[j,p])[],j=1..n)]; #Row(B,i).Column(A,p)
      end do
    end do;
    B:=C;
  end do;
  if not IsDirected(G) then # correct for double counting except for 2-cycles
        cycles:=map(i->i/2,cycles);
        if assigned(cycles[2]) then
            if cycles[2]=0 then # edgeless graphs
               evaln(cycles[2])
            else
               cycles[2]:=cycles[2]*2
            end if
        end if
  end if;
  eval(cycles);
end proc:

>

"G3:=CompleteGraph(5); cycles:=CycleCounts(G3): for i,cyc in eval(cycles) do i,cyc end do;add(eval(cycles)); "

>

memory used=138.66MiB, alloc change=0 bytes, cpu time=3.23s, real time=3.36s, gc time=156.25ms

>

Download ZeonsCycleCounts.mw

A linear algebra based version...

Answered: dharr 8445

June 21 2023

2 0

Smarter than Tarjan, perhaps, but slower:

>

Using the adjacency matrix and LinearAlgebra:-Generic to find cycles.

In this version, we use Arrays. In the modified adjacency matrix we just use and not .

>

Adj := AdjacencyMatrix(G); A := Matrix(n, n, proc (i, j) options operator, arrow; if Adj[i, j] = 1 then j else 0 end if end proc); B2 := Matrix(n, n, proc (i, j) options operator, arrow; if Adj[i, j] = 1 then {Array([i, j])} else {} end if end proc, datatype = set); B2, A

Matrix(%id = 36893489640244560108), Matrix(%id = 36893489640244559988)

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. WE ASSUME THE SECOND ONE HAS ONLY ONE VERTEX.

>

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.

>

Multiplication is adjoining the extra vertex. But if the new is the same as an existing one we have intersected ourself, but not necessarily made a simple cycle. So we want to stop adding vertices by making the matrix entry {}. We can identify it as a cycle if it intersects the first vertex.

>	map((x,y)->if y in x then NULL else (Array(x),=y) end if,a,b);

The unit for multiplication should leave the set of sets unchanged, so {Array([])} can be used. (but never called)

>	map((x,y)->if y in x then NULL else (Array(x),=y) end if,{Array([])},b); map((x,y)->if y in x then NULL else (Array(x),=y) end if,a,{Array([])}[]); #needs special case

And we should check that zero multiplied by is zero

>	map((x,y)->if y in x then NULL else (Array(x),=y) end if,{},b); map((x,y)->if y in x then NULL else (Array(x),=y) end if,a,{}); #treat as special case

Define these operations for the LinearAlgebraGeneric package:

>

cycles := table(); m := 0; F[`=`] := `=`; F[`/`] := `/`; F[`0`] := {}; F[`1`] := 1; F[`+`] := `union`

>

F[`*`]:=(x, y) -> if y = 0 or y={} or x = {} then {}
                    elif y = {Array([])} then x
                    else map(proc(u, w);
                               if w in u then
                                  if w = u[1] and w = min(u) then cycles[k][u] := NULL end if;                                                     return NULL
                               else
                                  return (Array(u),=w)
                               end if;
                             end proc, x, y)
                    end if:

Warning, (in anonymous procedure within F[*]) `cycles` is implicitly declared local

>

Matrix(%id = 36893489716538042956)

>

Matrix(%id = 36893489716513642604)

>

So here is a procedure for the cycles in a graph - each cycle appears twice.

>

Cycles:=proc(G::Graph)
  uses GraphTheory;
  local x,y,k,u,w,A,B,F,n,cycles;
  #if IsDirected(G) then error "graph must be undirected" end if;
  cycles:=table();
  F[`=`]:=`=`:F[`/`]:=`/`: # not required but have to be specified
  F[`0`]:={}:
  F[`1`]:={Array([])}: # never used
  F[`+`]:=`union`:     # can get arguments in any order
  F[`*`]:=(x, y)       # x from left Matrix, y from right Matrix
                 -> if y = 0 or x = {} then {}
                    #elif y = {Array([])} then x
                    else map(proc(u, w);
                               if w in u then
                                  if w = u[1] and w = min(u) then cycles[k][u] := NULL end if;
                                  return NULL
                               else
                                  return (Array(u),=w)
                               end if;
                             end proc, x, y)
                    end if;
  n:=NumberOfVertices(G);
  A:=Matrix(n,n, (i, j)-> if AdjacencyMatrix(G)[i, j] = 1 then j else 0 end if);
  B:= Matrix(n,n, (i, j)-> if AdjacencyMatrix(G)[i, j] = 1 then {Array([i,j])} else {} end if);
  for k from 2 to n do
    B:=LinearAlgebra:-Generic:-MatrixMatrixMultiply[F](B,A);
  end do;
  cycles:=map({indices},cycles,'nolist');
  #cycles[2]:=evaln(cycles[2]); # if don't include 2-cycles
  eval(cycles);
end proc:

>

memory used=0.57GiB, alloc change=-32.00MiB, cpu time=5.88s, real time=4.94s, gc time=1.89s

>

Download GenericCyclesArraysVerticesOnly2.mw

Tarjan's algorithm...

Answered: dharr 8445

June 21 2023

1 2

I looked into this a while ago, and found Tarjan's backtracking algorithm is quite efficient, implementing the stack with Maple's DEQueue (not using the double ended feature). I just implemented it as stated; it may be able to be improved (perhaps the stack clearing). This doesn't use the cycle basis though.

Tarjan's algorithm for finding cycles. R. Tarjan, Enumeration of the Elementary Circuits of a Directed Graph, SIAM J. Comput. 2 (1973) 211. doi: 10.1137/0202017

Procedure Tarjan (in startup code) accepts a Graph and outputs a table with all cycles of a graph, indexed by the size of the cycles.
For an undirected graph each cycle except the 2-cycles appears twice, once in each orientation (clockwise or counterclockwise).

Cycles of a given size are in sets, with each cycle gives as an ordered list of vertices in the cycle, with the vertices numbered 1,..n, in the same order as Vertices(G). The lowest numbered vertex is first in the cycle.

For counting cycles without producing the cycles themselves use procedure TarjanCounts, which produces a table indexed by the cycle sizes giving counts of cycles of that size. This time, each cycle counts once; there is no double counting for undirected graphs.

>

restart;

>	with(GraphTheory):

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

Cycles broken down by lengths.

>	cycles:=Tarjan(G):

>	for i,cycle in eval(cycles) do i, cycle; end do;

Total number of cycles

>	add(map(nops,cycles));

>	counts:=TarjanCounts(G): for i,cycle in eval(counts) do i, cycle; end do; add(eval(counts));

>

memory used=90.10MiB, alloc change=32.00MiB, cpu time=1.81s, real time=1.82s, gc time=125.00ms

Numbers with different lengths.

>

The 16 cycles of length 15

>

This version just counts cycles, and doesn't output them.

>

memory used=53.73MiB, alloc change=0 bytes, cpu time=812.00ms, real time=815.00ms, gc time=0ns

>

"for i,cycle in eval(counts) do i, cycle; end do; add(eval(counts)); "

>

Download Tarjan.mw

Code in startup code reproduced here (version that just counts is also there).

# Tarjan outputs a table with all cycles of a graph, indexed by the size of the cycles
# for an undirected graph each cycle except the 2-cycles appears twice, once for each orientation.
Tarjan:=proc(G::GRAPHLN);
  local V,AA,backtrack,cycles,marked_stack,point_stack,mark,A,u,i;
#  
  backtrack:=proc(s,v)
    local f,u,w,g;
    f:=false;
    push_back(point_stack,v);
    mark[v]:=true;
    push_back(marked_stack,v);
    for w in A[v] do
      if w<s then
        A[v]:=A[v] minus {w};
      elif w=s then
        cycles[numelems(point_stack)][convert(point_stack,list)]:=NULL; #record cycle
        f:=true;
      elif not mark[w] then
        g:=procname(s,w);
        f:=f or g;
      end if; 
    end do;
    if f then
      while back(marked_stack)<>v do
        u:=pop_back(marked_stack);
        mark[u]:=false
      end do;
      pop_back(marked_stack);
      mark[v]:=false;
    end if;
    pop_back(point_stack);
    return f;
  end proc:
#
  V:=GraphTheory:-NumberOfVertices(G):
  AA:=op(4,G):
  cycles:=table():
  marked_stack:=DEQueue():
  point_stack:=DEQueue():
  mark:=Vector(V,'fill'=false): 
  for i to V do
    A:=copy(AA);
    backtrack(i,i);
    while not empty(marked_stack) do
      u:=pop_back(marked_stack);
      mark[u]:=false;
    end do;
  end do;
  map({indices},cycles,'nolist');
end proc;

Need an accurate initial point to use hi...

Answered: dharr 8445

June 17 2023

2 3

>

restart;

>	fun:=(2(x+y+z)(sqrt(yz(z+x)(x+y))/(z+2x+y\ )+sqrt(zx(x+y)(y+z))/(x+2y+z)+sqrt(xy(y+z)(z+x))/(y\ +2z+x))-(9xyz/(x+y+z)+2(yz+zx+xy)))/(sqrt(xyz/(x+y+z))(x+y+z-sqrt(27xy*z/(x+y+z))));

(2*(x+y+z)*((y*z*(z+x)*(x+y))^(1/2)/(z+2*x+y)+(z*x*(x+y)*(y+z))^(1/2)/(x+2*y+z)+(x*y*(y+z)*(z+x))^(1/2)/(y+2*z+x))-9*x*y*z/(x+y+z)-2*x*y-2*z*x-2*y*z)/((x*y*z/(x+y+z))^(1/2)*(x+y+z-3*3^(1/2)*(x*y*z/(x+y+z))^(1/2)))

Apply a mix of the strategies of @mmcdara and @Axel Vogt
If we scale each variable by a, we do not change the value of fun.

>	simplify(eval(fun,{x=ax,y=ay,z=a*z})-fun) assuming a>0;

So we can arbitrarily set x=1, and reduce it to a to two-variable problem a la @mmcdara

>	f2:=eval(fun,x=1);

(2*(1+y+z)*((y*z*(z+1)*(1+y))^(1/2)/(z+2+y)+(z*(1+y)*(y+z))^(1/2)/(1+2*y+z)+(y*(y+z)*(z+1))^(1/2)/(y+2*z+1))-9*y*z/(1+y+z)-2*y-2*z-2*y*z)/((y*z/(1+y+z))^(1/2)*(1+y+z-3*3^(1/2)*(y*z/(1+y+z))^(1/2)))

As @Axel Vogt noted if we take z=1/y we achieve the alleged maximim at y=0 and y=infinity, suggesting we want y=0 and z=infinity or vice versa

>	f3:=eval(f2,z=1/y):plot(f3,y=0..infinity);

>	limit(f3,y=0,right);limit(f3,y=infinity);

So we take this as the initialpoint and see how well we will do

>	Optimization:-Maximize(f2, assume = nonnegative, initialpoint=[y=0,z=infinity]);

Error, (in Optimization:-NLPSolve) no improved point could be found

Well if we trust this, we have hit the jackpot. But let's try some values slightly off. The large z is unchanged

>	Optimization:-Maximize(f2, assume = nonnegative, initialpoint=[y=0.011,z=1e99]);

Any smaller y leads to

>	Optimization:-Maximize(f2, assume = nonnegative, initialpoint=[y=0.010,z=1e99]);

Error, (in Optimization:-NLPSolve) no improved point could be found

This calculation is apparently done with the NAG routines (From the package overview help: "The package takes advantage of built-in library routines provided by the Numerical Algorithms Group (NAG)"), which I took to mean it was done at hardware precision. Therefore I was expecting it to not to improve once Digits was higher than the evalhf(Digits) value (here 15). But it does seem possible now we are close to the right answer:

>	interface(warnlevel=0): #suppress "Warning, undefined value encountered" Digits:=200:Optimization:-Maximize(f2, assume = nonnegative, initialpoint=[y=0.011,z=1e99]);

>

Download maximize.mw

not all x,y,z,w give solutions...

Answered: dharr 8445

June 15 2023

1 0

>

eq := a+b+c = n*x; eq2 := a+3*b+c = n*y; eq3 := a+b+5*c = n*z; eq4 := 4*a+8*b = n*w

>

Matrix(%id = 36893489914357921060), Vector[column](%id = 36893489914357920940)

Consider the case used by @tomleslie

>

Matrix(%id = 36893489914357910812)

Since the Rank is 3, there is an infinite set of one-parameter solutions, as @tomleslie found from isolve, and they are integer

>

Vector[column](%id = 36893489914357963100)

If the rank is 4 then the only solution will be a=b=c=n=0., which doesn't have n positive.

So the only sets of {w,x,y,z} that have solutions (integer or otherwise) will be those that make the last column a linear combination of the other 3 columns. The solutions are then the coefficients of the linear combination. For example we can add the columns

>

Matrix(%id = 36893489914345976044)

and so the solution is (by construction) {a=b=c=n=1} or any multiple of this.

>

Vector[column](%id = 36893489914345981340)

>

Download LinearSystem.mw

one way...

Answered: dharr 8445

June 12 2023

1 0

Here's one way.

Download cis.mw

Heaviside or piecewise...

Answered: dharr 8445

June 07 2023

1 1

Here are two suggestions.

Download Heaviside.mw

slightly different...

Answered: dharr 8445

June 03 2023

0 2

For me (2023.0, Windows 10), the second one hangs, but if I use the stop icon (after 250 s) it quickly ends and prints "done", and then I can continue to use Maple as normally.

Using dsolve with NLPSolve...

Answered: dharr 8445

May 31 2023

1 7

Here's how I do least squares fits with numerical solutions from dsolve. I bypass the parameters facility of dsolve. I've only tested this on the derivative-free nonlinearsimplex method [but see below]. I use the matrix method call to NLPSolve; for the procedure method it seemed I had to hard-code the number of parameters in the internal procedure SS.

Edits: it seems to work with sqp, but the initial point had to be closer. Updated with streamlined version.

>

restart;

Test DE with parameters in DE and initial conditions

>	de:={diff(y(t),t,t)=-a^2y(t),D(y)(0)=ba,y(0)=c};

${diff(diff(y(t), t), t) = -a^2*y(t), y(0) = c, (D(y))(0) = b*a}$

Analytical answer

>	ans:=rhs(dsolve(de,y(t)));

Generate signal for some test parameters

>	params:={a=1.1,b=2.,c=3.}; sig:=unapply(eval(ans,params),t); tt:=Vector([seq(t,t=0..6,0.1)]): yy:=map(sig,tt):

>	#plot(tt,yy,style=point);

Routine to fit to differential equations with parameters and one measured concentration (first one given in vars).

Outputs the residual sum of the squares and the parameter values found in the same order as in the initial guesses.
Uses NLPSolve's matrix input option

>

defit:=proc(tims::~Vector,expt1::~Vector,des::set,vars::list,inivarsarg::list)
  local inivals,gparams,RSS;
  gparams:=map(lhs,inivarsarg);
  inivals:=map(rhs,inivarsarg);
  RSS:=proc(V)
    local sol;
    if not type(V,'Vector'(numeric)) then return 'procname'(args) end if;
    sol:=dsolve(eval(des,gparams=~convert(V,list)),vars,numeric,output=Array(tims),method=rosenbrock);
    add((expt1-sol[2,1][..,2])^~2);
  end proc:
  Optimization:-NLPSolve(nops(inivals),RSS,initialpoint=Vector(inivals),method=sqp);
end proc:

>	defit(tt,yy,de,[y(t)],[a=1.5,b=2.5,c=3.3]);

$[9.12450471859734705*10^(-13), Vector(3, {(1) = 1.1000000620919381, (2) = 1.9999999531931567, (3) = 2.9999998723903567})]$

>

Download defit.mw

assign...

Answered: dharr 8445

May 27 2023

0 0

Download assign.mw

You can do the assign as you read a line by putting a line such as "splog=0.56" into a variable str, and then use

assign(parse(str));

convex hull...

Answered: dharr 8445

May 21 2023

5 1

In this case, the points of the polyhedron and its convex hull are the same, so ConvexHull in the ComputationalGeometry package works

Download volume.mw

x range, numpoints...

Answered: dharr 8445

May 20 2023

1 0

If you work out the exact x range for the spacecurve it works, e,g, for the circle it is just x=-3..3. For the ellipse see the attached. Alternatively, you can just use numpoints=10000

SectionsConiquesTest.mw.

different T...

Answered: dharr 8445

May 17 2023

2 8

The T you used in ConvertIn(T) is not the same as the internal one. The easiest way to get around this is always to provide an irreducible polyomial to GF, then the variable you use there can be used.

Download GF.mw

CUDA...

Answered: dharr 8445

May 16 2023

3 30

CUDA only helps speed up numerical floating point matrix multilpication, which will not help your symbolic solution of polynomial equations. See ?CUDA,supported_routines