====== TQmat ======

Notes:
  * the network matrices are considered two-mode

<code>
# operations with temporal quantities (tau = 0)
# March 10, 2014

import pickle

a = [(1,5,2), (6,8,1), (11,12,3), (14,16,2), (17,18,5), (19,20,1)]
b = [(2,3,4), (4,7,3), (9,10,2), (13,15,5), (16,21,1)]
inf = float("inf"); one = 1; zero = 0
N = []
E = [(1,inf,one)]   # the third component is the unit in semiring

def TQget(S):
   try: return(next(S))
   except StopIteration: return((inf,inf,0))

def plus(a,b): return(a+b) 
def times(a,b): return(a*b) 

def TQstandard(a):
   if len(a) == 0: return(a)
   fc = inf; c = []
   for (sa,fa,va) in a:
      if fc == sa:
         if vc == va: fc = fa
         else:
            if sc != fc: c.append((sc,fc,vc))
            vc = va; sc = sa; fc = fa
      else:
         if (fc != inf) and (sc != fc): c.append((sc,fc,vc))
         sc = sa; fc = fa; vc = va
   if sc != fc: c.append((sc,fc,vc))
   return(c)

def TQbinary(a):
   c = []
   for (sa,fa,va) in a:
      if va > 0: c.append((sa,fa,1))
   return(TQstandard(c))

def TQinvert(a):
   c = []
   for (sa,fa,va) in a:
      c.append((sa,fa,1/va))
   return(c)

def TQcentral(v,C):
   cc = c = []; n = len(C)
   for t in range(n):
      if t!=v:
         c = TQsum(c,TQinvert(C[v][t]))         
   for (s,f,v) in c: cc.append((s,f,v/(n-1)))
   return(cc)
                               
def TQsum(a,b):
   if len(a) == 0: return(b)
   if len(b) == 0: return(a)
   c = []; A = a.__iter__(); B = b.__iter__()
   (sa,fa,va) = TQget(A); (sb,fb,vb) = TQget(B)
   while (sa<inf) or (sb<inf):
      if sa < sb:
         sc = sa; vc = va
         if sb < fa: fc = sb; sa = sb
         else: fc = fa; (sa,fa,va) = TQget(A)        
         c.append((sc,fc,vc))
      elif sa == sb:
         sc = sa; fc = min(fa,fb); vc = plus(va,vb)
         c.append((sc,fc,vc))
         sa = sb = fc; fA = fa
         if fA <= fb: (sa,fa,va) = TQget(A)
         if fb <= fA: (sb,fb,vb) = TQget(B)
      else:
         sc = sb; vc = vb
         if sa < fb: fc = sa; sb = sa
         else: fc = fb; (sb,fb,vb) = TQget(B)        
         c.append((sc,fc,vc))
   return(TQstandard(c))

def TQprod(a,b):
   if len(a)*len(b) == 0: return([])
   c = []; A = a.__iter__(); B = b.__iter__()
   (sa,fa,va) = TQget(A); (sb,fb,vb) = TQget(B)
   while (sa<inf) or (sb<inf):
      if fa <= sb: (sa,fa,va) = TQget(A)
      elif fb <= sa: (sb,fb,vb) = TQget(B)
      else:
         sc = max(sa,sb); fc = min(fa,fb); vc = times(va,vb)
         c.append((sc,fc,vc))
         if fc == fa: (sa,fa,va) = TQget(A)
         if fc == fb: (sb,fb,vb) = TQget(B)          
   return(TQstandard(c))         

def TQlist(a): # ???
   for (sa,fa,va) in a:
      c.append((sa,fa,1/va))
   return(c)

def VecOne(n):
   return([ E for i in range(n) ])

def MatList(M,names=None,all=True):
   n = len(M)
   for i in range(n):
      for j in range(n):
         if len(M[i][j])>0:
            if names==None: print('(',i+1,',',j+1,') =',M[i][j])
            else:
               print('(',names[i],',',names[j],') =\n  ',M[i][j])

def MatClosure(R):
   n = len(R); C = R
   for k in range(n):
      for u in range(n):
         for v in range(n): 
            C[u][v] = TQsum(C[u][v], TQprod(C[u][k],C[k][v]))
      C[k][k] = TQsum(E,C[k][k])
   return(C)

def MatBin(A):
   nr = len(A); nc = len(A[0])
   B = [[ [] for i in range(nc)] for j in range(nr)]
   for u in range(nr):
      for v in range(nc):
         B[u][v] = TQbinary(A[u][v])
   return(B) 

def sumSyM(C,Rows,Cols,diag=False):
# extend for Rows and Cols are temporal partitions
   n = len(C); s = []
   for u in Rows:
      for v in Cols:
         if u>v: s = TQsum(s,C[v][u]) 
         elif u<v: s = TQsum(s,C[u][v]) 
         elif diag: s = TQsum(s,C[v][u])
   return(s)

def MatConst(A,const):
   nr = len(A); nc = len(A[0])
   B = [[ [] for i in range(nc)] for j in range(nr)]
   for u in range(nr):
      for v in range(nc):
         B[u][v] = TQprod(const,A[u][v])
   return(B) 

def MatSum(A,B):
   nra = len(A); nca = len(A[0]); nrb = len(B); ncb = len(B[0]);
   if nra!=nrb: return(None)
   if nca!=ncb: return(None)
   C = [[ [] for i in range(nca)] for j in range(nra)]
   for u in range(nra):
      for v in range(nca):
         C[u][v] = TQsum(A[u][v],B[u][v])
   return(C)

def MatVecLeft(A,x): 
   nr = len(A); nc = len(A[0]); nv = len(x)
   if nv!=nr: return(None) 
   y = [ [] for i in range(nc)] 
   for v in range(nc):
      s = [] 
      for u in range(nr):
         s = TQsum(s,TQprod(x[u],A[u][v]))
      y[v] = s
   return(y)

def MatVecRight(A,x): 
   nr = len(A); nc = len(A[0]); nv = len(x)
   if nv!=nc: return(None) 
   y = [ [] for i in range(nr)] 
   for u in range(nr):
      s = [] 
      for v in range(nc):
         s = TQsum(s,TQprod(A[u][v],x[v]))
      y[u] = s
   return(y)

def MatVecRow(A,x): 
   nr = len(A); nc = len(A[0]); nv = len(x)
   if nv!=nr: return(None) 
   B = [[ [] for i in range(nc)] for j in range(nr)]
   for u in range(nr):
      for v in range(nc):
         B[u][v] = TQprod(x[u],A[u][v])
   return(B)

def MatVecCol(A,x): 
   nr = len(A); nc = len(A[0]); nv = len(x)
   if nv!=nr: return(None) 
   B = [[ [] for i in range(nc)] for j in range(nr)]
   for u in range(nr):
      for v in range(nc):
         B[u][v] = TQprod(A[u][v],x[v])
   return(B) 

def MatProd(A,B):
   nra = len(A); nca = len(A[0]); nrb = len(B); ncb = len(B[0]);
   if nca!=nrb: return(None)
   C = [[ [] for i in range(ncb)] for j in range(nra)]
   for u in range(nra):
      for v in range(ncb):
         s = [] 
         for t in range(nca):
            s = TQsum(s,TQprod(A[u][t],B[t][v]))
         C[u][v] = s
   return(C)

def MatProdDiag(A,B):
   nra = len(A); nca = len(A[0]); nrb = len(B); ncb = len(B[0]);
   if nca!=nrb: return(None)
   if nra!=ncb: return(None)
   x = [ [] for u in range(nra)] 
   for u in range(nra):
      s = [] 
      for v in range(nca):
         s = TQsum(s,TQprod(A[u][v],B[v][u]))
      x[u] = s
   return(x)

def MatSetDiag(A,const):
   nr = len(A); nc = len(A[0])    
   if nr!=nc: return(None)
   B = A
   for u in range(nr): B[u][u] = const
   return(B)

def MatSym(A):
   nr = len(A); nc = len(A[0])
   if nr!=nc: return(None)
   B = A
   for u in range(nr-1):
      for v in range(u+1,nr):
         B[u][v] = TQsum(A[u][v],A[v][u])
         B[v][u] = B[u][v]
   return(B)

def inDeg(A):
   return(MatVecLeft(MatBin(A),VecOne(len(A)))) 
   
def outDeg(A):
   return(MatVecRight(MatBin(A),VecOne(len(A[0])))) 
   
def clusCoef(A,type=1):
# type = 1 - standard clustering coefficient
# type = 2 - corrected clustering coefficient / temporal degMax 
# type = 3 - corrected clustering coefficient / overall degMax
# 9-10. april 2014
   nr = len(A); nc = len(A[0])
   if nr!=nc: return(None)
   Ab = MatSetDiag(MatBin(A),N)
   S = MatSym(Ab)
   deg = MatVecRight(MatBin(S),VecOne(nr))
   if type == 1:
      inv = []
      for d in deg:
         e = []
         for (s,f,v) in d:
            if v > 1: e.append((s,f,1/v/(v-1)))
         else: e.append((s,f,0))
         inv.append(e)
#      print("Invert")
#      for (i,e) in enumerate(inv): print(i+1,e)
      tri = MatProdDiag(MatProd(S,Ab),S)
      cc = []
      for (i,t) in zip(inv,tri): cc.append(TQprod(i,t))
      return(cc)
   else: return(None)
</code>