> restart:with(combinat):with(linalg):

Warning, the protected name Chi has been redefined and unprotected

Warning, the name fibonacci has been rebound

Warning, the protected names norm and trace have been redefined and unprotected

The procedure matroid tests to see whether a given r by r+n matrix M represents a matroid of rank r.

> matroid:=proc(r,n,M) local A,AA,i,j,k,l,m,z,v,S,test,bad,floop,T,V,W,dun,dup,sum1,sum2,modsum,ind,findind,total,prob1,prob2:

> global mst,numtests,flagmat,lastspot:#numtests:=numtests+1:#flagvec:=vector(n,0):

The variables r (rank), n (number of chords), and M (the matrix representing the matroid) are passed to this procedure.

The local variables are as follows:

"test" is a vector containing the lengths of the circuits computed so far.  A is the set of lengths of circuits.  (The components of test are the members of A.)

AA is not used; it was diagnostic.  i, j, k, l, m, and z are dummy variables for summation.  v is a vector of column vectors (a matrix); the columns are the columns of M to the right of the r x r identity matrix.  S is the set of these vectors.

"bad" is an n by (n choose [n/2]) matrix.  It is used to record which subsets of columns do not yield circuits.  

"floop" is a dummy looping variable for flagging

T is a list of sets.  Each of these sets has all subsets of S of a certain size, starting with 2.  The sets are indexed by the size of subsets they contain.  For example, T[5] is the list of all 5-element subsets of S.  T[i][j] is the jth subset in the list of all i-element subsets. It contains vectors from S, which are columns of M to the right of the identity.  We will need to examine these to decide whether they, along with some columns of the identity portion of M, form a circuit.

V is a matrix of lists, more or less.  V looks at the T[i][j], which are sets of columns, and finds subsets of those.  V[i,j,k] is the list of all k-element subsets of columns drawn from a particular subset of columns T[i,j].  In order to determine whether the columns appearing in T[i,j] (along with appropriate columns from I_r) form a circuit, we need to ensure that no proper subset of these columns is dependent since circuits must be minimally dependent.  This is the purpose of V[i,j,k].

V[i,j,k,l] and W[i,j,k,l] are lists.  V[i,j,k,l] is the l-th (ell-th) member of the list of the columns in one of the k-element subsets of T[i][j].  W[i,j,k,l] is a list of the columns not included in V[i,j,k,l]; it is the complement of V[i,j,k,l] in T[i][j].

dup is not used in this incarnation.

MAPLE will look at every set and its complement.  However, since every set is the complement of its complement, I need a way to avoid looking at the same set twice.  dun keeps track of which sets have already been examined and exits a loop when it spots a duplicate.

sum1 is the sum of the vectors in V[i,j,k,l] and sum2 is the sum of the vectors in W[i,j,k,l].  modsum is the sum of these two vectors mod 2.

ind is an index, used exclusively in test. findind is an index used to compare previously computed circuit sized with the current circuit size.  When a circuit size is duplicated, the offending f-vectors are recorded in prob1 and prob2.

total is the sum of the vectors in T[i][j] mod 2.  If the procedure gets to this point, then total is actually used to compute the length of the circuit: it is the number of 1's in total plus i (the number of vectors).  This number is then stored in test.

mst is diagnostic.  numtests keeps track of how many times the procedure has been applied; it is disabled for this one-at-a-time version.

> mst:=time():flagmat:=matrix(n,r-1,0):

> ind:=0:S:={}:A:={}:dup:={}:AA:=vector(2^n,0):

> #for i from 1 to r+1 do dup:=dup union{col(M,i)}:od:

The loop below extracts the columns of M to the right of the r x r identity matrix; v[i] is the ith such column.

>

> for i from 1 to n do

> v[i]:=col(M,r+i):#if v[i] in dup or v[i] in S or v[i]=vector(r,0) then return:fi:

> S:=S union {v[i]}:

> od:

>

> test:=vector(r-1,0):

>

testing and measuring subsets

> bad:=matrix(n,binomial(n,floor(n/2)),0):

>

>

> for j from 1 to n do

> ind:=ind+1:if ind >= r then return: fi:test[ind]:=dotprod(v[1],v[j])+1:flagmat[j,ind]:=1:

> if test[ind] in A then flagvec[j]:=1:return:fi:

> A:=A union {test[ind]}:

> if 1 in A or 2 in A then return:fi:

> od:

>

> for i from 2 to n do

>

> T[i]:=convert(choose(S,i),list):#print(T[i]);

>

> for j from 1 to binomial(n,i) do

>

>

I need to now analyze each subset for compatibility.  j becomes the leading index.

> dun:={}:

> for k from 1 to floor(i/2) do

>

> V[i,j,k]:=convert(choose(T[i][j],k),list):

>

Vijk is the list of k element subsets of Tij.  I need to compare each of these to its complement in Tij.

> for l from 1 to binomial(i,k) do

> W[i,j,k,l]:=convert(T[i][j] minus V[i,j,k][l],list):

> V[i,j,k,l]:=convert(V[i,j,k][l],list):

>

Wijk is the set of complements of the Vijk in Tij.

> if V[i,j,k,l] in dun or W[i,j,k,l] in dun then next:next:next:fi:

> dun:=dun union {V[i,j,k,l],W[i,j,k,l]}:

>

> sum1:=vector(r,0):sum2:=vector(r,0):

> for m from 1 to k do

> sum1:=matadd(sum1,V[i,j,k,l][m]):

> od:

>

> for m from 1 to i-k do

> sum2:=matadd(sum2,W[i,j,k,l][m]):

> od:

>

>

>

This tests for incompatible subsets.

> for z from 1 to r do sum1[z]:=irem(sum1[z],2): sum2[z]:=irem(sum2[z],2): od:

> modsum:=matadd(sum1,sum2):

> for z from 1 to r do modsum[z]:=irem(modsum[z],2):od:

> if dotprod(v[1],sum1)+dotprod(v[1],sum2)=dotprod(v[1],modsum) then bad[i,j]:=1:next:next:next:

The above line does the following: The f_k in T[i][j] together with the appropriate e_i (chosen minimally) are dependent.  If some subset of them, along with an appropriate subset of the chosen e_i, are also dependent, then T[i][j] does not correspond to a circuit since it is not *minimaly* dependent.  If the sum of these f_k has no 1's in common with the sum of the other f_k, this will happen -- there is no interference from the other f_k.  This portion of the procedure determines this by adding them without reducing mod 2 and comparing to the mod 2 sum.  If there is overlap, these will not be equal because the former will contain a 2.  If there is not overlap, these two will be equal and this set will be thrown out.  Otherwise, we have a circuit and its length is computed below.

> fi:

>

> od:

> od:

>

>

The following computes the length of the circuit (if it is one -- indicated by bad[i,j]=0).

>

> if bad[i,j]=0 then total:=vector(r,0):

> for m from 1 to i do

> total:=matadd(total,T[i][j][m]):

> od:

> for z from 1 to r do

> total[z]:=irem(total[z],2):

> od:

> ind:=ind+1:if ind=r then return;fi;#print(r, `is too low a rank`); return;fi;

> test[ind]:=dotprod(v[1],total)+i:#T[i][j]:

>

>

>

> for floop from 1 to n do

> if v[floop] in T[i][j] then flagmat[floop,ind]:=1:fi:od:

>

>

> if test[ind] in A then

> for findind from 1 to ind-1 do

> if test[findind]=test[ind] then prob1:=col(flagmat,findind):prob2:=col(flagmat,ind):fi:

> od:

> for floop from 1 to n do

> if prob1[floop]=1 or prob2[floop]=1 then lastspot:=floop:fi:

> od:

>

> #print(T[i][j],lastspot,flagmat,test,prob1,prob2,M);

>

> return:fi:

>

>

> if 1 in A or 2 in A then return:fi:

> A:=A union {test[ind]}:#flagvec[ind]:=test[ind]:

> fi:

>

>

> od:od:if test[r-1]=0 then return;fi;#print(r, `is too high a rank`); return;fi;

>

>

> print(r,test,M);

> end proc:

>

Warning, `flagvec` is implicitly declared local to procedure `matroid`

> M:=matrix(24,29,0):

> for i from 1 to 24 do

> M[i,25]:=1:M[i,i]:=1:

> od:

> M[1,26]:=1:M[2,26]:=1:M[3,27]:=1:M[4,27]:=1:M[5,27]:=1:M[1,28]:=1:

> for i from 7 to 13 do

> M[i,28]:=1:M[i-5,29]:=1:

> od:

> M[8,29]:=0:

> print(M);

> matroid(24,5,M);

>