lcz

620 Reputation

10 Badges

3 years, 89 days
changsha, China

MaplePrimes Activity


These are answers submitted by lcz

My advice is to look at the source codes of MathChem and write them yourself. I saw the introduction of package MathChem you mentioned.  MathChem is a python package for calculating topological Indices. This is actually a similar question I've asked before, but I talked about the igraph package.

Another similar discussion is as below.

I tried to install this package in python3, but unfortunately the run always failed. So I looked at the source codes, and it seemed to be based on Python 2.  I looked at  some functions in MathChem , and most of them could have been rewritten in Maple or made use of the Maple built-in functions.  

For example,  the spectral radius of adjacency matrix of a graph is what we talked about a while ago.  Carl Love,  dharr and acer  have an equally brilliant analysis of the problem.

SR := proc(G::Graph)
  local M := GraphTheory:-AdjacencyMatrix(G);
  max(abs~(LinearAlgebra:-Eigenvalues(
    rtable(`if`([rtable_indfns(M)]=[':-symmetric'],
           ':-symmetric',NULL), M,
           'datatype'=':-hfloat','storage'=':-rectangular',
           'order'=':-Fortran_order'))));
end proc:

  See more details in https://www.mapleprimes.com/questions/233883-How-To-Quickly-Find--The-Spectral-Radius-Of-A-Matrix-.

Another example, we can calculate the first Zagreb index and  second Zagreb index of a graph by maple.

with(GraphTheory):with(SpecialGraphs):
G:=PetersenGraph()
Zagrebindex_1:=g-> add((GraphTheory:-DegreeSequence(g))^~(2));
Zagrebindex_2:=proc(g)
local d,i;
d:=Edges(g);
add([seq(Degree(g,d[i][1])*Degree(g,d[i][2]), i=1..NumberOfEdges(g))])
end proc:
Zagrebindex_1(G),Zagrebindex_2(G)

90, 135

In MathChem, their codes are written in a similar way.

def zagreb_m1_index(self):
        """ Zagreb M1 Index """    
        return sum(map(lambda d: d**2, self.degrees()))
def zagreb_m2_index(self):
        """ Zagreb M2 Index 
        
        The molecular graph must contain at least one edge, otherwise the function Return False
        Zagreb M2 Index is a special case of Connectivity Index with power = 1"""
        return sum( map(lambda (e1, e2): self.degrees()[e1]*self.degrees()[e2] , self.edges()) )            
    
   

During my master's period, I once did research on chemical graph theory. Generally speaking, it is not too difficult to write codes for chemical index of graph. The difficulty is  calculating some parameters of graphs themselves.  For example, whether a graph is Hamiltonian, or to calculate  the crossing number  of a graph, and so on.  (They are NP hard. )  Even for linear time algorithms, such as testing whether a graph is planar, it is difficult to write code without specialized knowledge, but fortunately it is often readily available.

The unweighted  version was previously provided in https://www.mapleprimes.com/questions/233695-Solve-The-Resistor-Grid-Of-1-Ohm-. Here, for the first question, we provide a  weighted graph version. Relevant theories can be referred to this paper.

  • Bapat R B. Resistance matrix of a weighted graph[J]. MATCH Commun. Math. Comput. Chem, 2004, 50(02).

As for whether the author discovered it for the first time, I haven't checked it carefully.

We just need to change the  laplacian matrix of  unweighted graph to the laplacian matrix of the weighted graph. So we improved the codes from last time as follows.

Laplaceweightmatrix:=proc(g::Graph)
local M,M1,M2;
M:=GraphTheory:-WeightMatrix(g):
M1:=Matrix(upperbound(M), (i,j)-> `if`(M(i,j)<>0, -1/M(i,j), 0)):
M2:=M1+LinearAlgebra:-DiagonalMatrix(-MTM:-sum(M1,1));
end proc:
ResistanceMatrix:= (g::Graph)-> 
    (M-> Matrix(
            upperbound(M), (i,j)-> M[i,i]+M[j,j], 'shape' = 'symmetric'
         ) - 2*M
    )(LinearAlgebra:-MatrixInverse(Laplaceweightmatrix(g)
          , 'method'= 'pseudo'
     ))
:
with(GraphTheory):
G:=Graph({[{"000","400"},4],[{"400","450"},5],[{"450","050"},4],[{"050","000"},5],
          [{"003","403"},4],[{"403","453"},5],[{"453","053"},4],[{"053","003"},5],
          [{"000","003"},3],[{"400","403"},3],[{"450","453"},3],[{"050","053"},3]}):
S:=ResistanceMatrix(G):
S[1,8] #156/47

 

PS: Here I have not studied the corresponding relationship between graph vertex index and adjacency matrix index thoroughly.

However,  it seems that Rows and columns of the adjacency matrix follow the order given by Vertices:

Vertices(G);

["000", "003", "050", "053", "400", "403", "450", "453"]

In graph theory, the resistance distance between two vertices of a simple connected graph, G, is equal to the resistance between two equivalent points on an electrical network, constructed so as to correspond to G, with each edge being replaced by a 1 ohm resistance

SEE https://en.wikipedia.org/wiki/Resistance_distance

PS: I remember that even though each edge is not set to a unit resistance, there is a corresponding theoretical calculation formula. 

So we can write the following function:

ResistanceMatrix:=proc(g::Graph)
   local L,ML,n, i, j, R;
     n:=GraphTheory:-NumberOfVertices(g);
     L:=GraphTheory:-LaplacianMatrix(g):
     ML:=LinearAlgebra:-MatrixInverse(L, method = pseudo):
         for i from 1 to n do
            for j from i to n do
               R:=Matrix(n,n,(i,j)->ML[i,i]+ML[j,j]-2*ML[i,j],shape = symmetric);
            od:
         od:
     R;
end proc:
with(SpecialGraphs):
S := SoccerBallGraph():
R:=ResistanceMatrix(S):
R[1,31];
R[1,2];

17/11;
16273/25080

Note that finding the  pseudo inverse of the matrix is ​​not efficient, so this method seems to be limited to cases with fewer vertices or only has theoretical value. Is there a faster solution? I don't know yet.

See some talks in https://www.mapleprimes.com/maplesoftblog/208427-Google-Maps-And-Geocoding-For-Maple 

First you have to install this extra package, but even so, it looks like this function is not easy to use now.

Page 1 of 1