====== NetsBM - blockmodeling based on Nets ======

<code python>

wdir = "C:/Users/batagelj/work/projects/BM/py"
gdir = "c:/users/batagelj/work/python/graph/Nets"
ddir = "C:/Users/batagelj/work/projects/BM/data"

import sys, os, re, json
sys.path = [gdir]+sys.path; os.chdir(wdir)
from Nets import Network as N
import numpy as np
from copy import copy, deepcopy
from timeit import default_timer as timer
from datetime import datetime

def P(net,C,Types):
    n = net._info['nNodes']; nC = Types.shape[0]; size = np.zeros(nC,dtype=int)
    for i in C: size[i] += 1
    dCnt = np.zeros(nC); Cnt = np.zeros((nC,nC))
  # Sum = np.zeros((nC,nC)); rSum = np.zeros((nC,n)); cSum = np.zeros((nC,n))
    rCnt = np.zeros((nC,n)); cCnt = np.zeros((nC,n))
    rPos = np.zeros((nC,nC)); cPos = np.zeros((nC,nC))
    error = np.zeros((nC,nC)); btype = np.zeros((nC,nC),dtype=int)
    for a in net.links():
        u,v,d,r,W = net.link(a); w = W['w']; p = u-1; q = v-1
        Cnt[C[p],C[q]] += 1
        rCnt[C[q],p] += 1; cCnt[C[p],q] += 1
        if p==q: dCnt[C[p]] += 1
    rB = (rCnt>0).astype(int); cB = (cCnt>0).astype(int)        
    for i in range(nC):
        for s in range(n):
            rPos[i,C[s]] += rB[i,s]; cPos[C[s],i] += cB[i,s]
    for i in range(nC):
        for j in range(nC):
            err = np.inf; typ = XXX
            if NUL in Types[i,j]: 
                ert = Cnt[i,j]
                if i==j: ert += min(0,size[i]-2*dCnt[i])
                if ert < err: err = ert; typ = NUL
            if COM in Types[i,j]: 
                ert = size[i]*size[j] - Cnt[i,j]
                if i==j: ert += min(0,2*dCnt[i]-size[i])
                if ert < err: err = ert; typ = COM
            if REG in Types[i,j]: 
                ert = (size[j]-cPos[i,j])*size[i] +\
                      (size[i]-rPos[i,j])*cPos[i,j]
                if ert < err: err = ert; typ = REG
            if RRE in Types[i,j]: 
                ert = (size[i]-rPos[i,j])*size[j]
                if ert < err: err = ert; typ = RRE
            if CRE in Types[i,j]: 
                ert = (size[j]-cPos[i,j])*size[i] 
                if ert < err: err = ert; typ = CRE
            error[i,j] = err; btype[i,j] = typ
    return (np.sum(error), error, btype)

def optimize(net,C,Types,minSize):
    nC = Types.shape[0]; size = np.zeros(nC,dtype=int)
    for i in C: size[i] += 1
    if np.min(size)<minSize: return(C,np.inf,None,None)
    found = True; Clu = np.array(C,dtype=int); Cm = Clu.copy()
    Err, Qm, Tm = P(net,Clu,Types)
    while found:
        found = False
        for a in net.arcs():
            u,v,d,r,W = net.link(a); w = W['w']
            p = u-1; q = v-1; Cp = Clu[p]; Cq = Clu[q]
            if not(Cp==Cq):
                Clu[p] = Cq; size[Cp] += -1; size[Cq] += 1
                if size[Cp]<minSize: err = np.inf
                else: err, Q, T = P(net,Clu,Types)
                if err<Err:
                    Cm = Clu.copy(); Err = err; found = True
                    Qm = Q.copy(); Tm = T.copy()
                else:
                    Clu[p] = Cp; Clu[q] = Cp; size[Cp] += 2; size[Cq] += -2
                    if size[Cq]<minSize: err = np.inf
                    else: err, Q, T = P(net,Clu,Types)
                    if err<Err:
                        Cm = Clu.copy(); Err = err; found = True
                        Qm = Q.copy(); Tm = T.copy()
                    else:
                        Clu[p] = Cq; size[Cp] += -1; size[Cq] += 1
                        err, Q, T = P(net,Clu,Types)
                        if err<Err:
                            Cm = Clu.copy(); Err = err; found = True
                            Qm = Q.copy(); Tm = T.copy()
                        else: Clu[p] = Cp; Clu[q] = Cq 
    return (Cm, Err, Qm, Tm)

def searchBM(nSteps=10,Pbest=np.inf,minSize=1):
    print('netsBM:', datetime.now())
    start = timer(); n = net._info['nNodes']; found = False; bad = 0
    iD = np.array([net.inDegree(u) for u in net.nodes()],dtype='I')        
    oD = np.array([net.outDegree(u) for u in net.nodes()],dtype='I')
    ind = np.logical_and(iD<2,oD<2); nSpec = np.sum(ind); spec = np.any(ind)
    for k in range(nSteps):
        j = k % 3
        if j==0: C = (nC*rng.random(n)).astype(int)
        elif j==1: C = C0.copy(); rng.shuffle(C)
        elif j==2:
            C = spec+((nC-spec)*rng.random(n)).astype(int)
            if spec: C[ind] = 0
        R = optimize(net,C,Types,minSize); Er = R[1]
        if Er == np.inf: bad += 1
        elif Er < Pbest:
            Pbest = Er; Rbest = deepcopy(R); found = True
            if n > 50: print(k+1,j,Er)
            else: print(k+1,j,Er,"\n  C :",C,"\n  C*:",R[0])
    end = timer(); print("time:",end - start,"s")
    if bad>0: print(bad,"bad initial partitions")
    if found: return Rbest

ddir = "C:/Users/batagelj/work/Python/graph/Nets/BM/dat"
net = N.loadPajek(ddir+"/class.net"); net.Info()

XXX = 0; NUL = 1; COM = 2; REG = 3; RRE = 4; CRE = 5
rng = np.random.default_rng()
n = net._info['nNodes']; nC = 3
C0 = np.array([0]*3 + [1]*3 + [2]*9, dtype='I')
</code>

===== Regular BM =====

<code python>
========== RESTART: C:/Users/batagelj/work/projects/BM/py/netsBM.py ==========
NODES:  n = 15  mode = 1
ARCS:  multirel = False

network:  test 
 Test 
simple= False  directed= True  org= 1  mode= 1  multirel= False  temporal= False 
nodes= 15  links= 56  arcs= 56  edges= 0
>>> Types = np.array([{NUL,COM,REG}]*(nC*nC),dtype='O').reshape(nC,nC)
>>> Reg = searchBM(500,minSize=2)
netsBM: 2020-08-15 00:09:39.810269
1 0 20.0 
  C : [2 2 0 1 0 1 1 2 2 1 2 2 1 2 2] 
  C*: [2 1 0 1 1 1 1 1 0 1 2 1 1 2 1]
3 2 19.0 
  C : [0 1 0 2 2 0 2 2 2 2 0 1 0 2 1] 
  C*: [0 1 2 2 2 2 1 2 2 2 0 1 2 1 1]
6 2 17.0 
  C : [1 2 2 2 2 0 1 1 0 0 1 0 2 1 2] 
  C*: [1 2 2 2 2 2 2 2 0 0 1 0 2 1 2]
11 1 15.0 
  C : [1 0 0 2 2 2 2 1 0 2 2 2 2 1 2] 
  C*: [1 2 2 2 2 2 2 2 2 2 0 2 2 1 0]
22 0 12.0 
  C : [0 2 0 1 2 1 0 2 1 2 2 0 1 2 0] 
  C*: [2 0 1 1 1 1 0 1 1 1 2 0 1 2 0]
97 0 11.0 
  C : [0 0 2 0 1 2 1 2 0 2 0 1 1 0 1] 
  C*: [0 1 2 2 1 2 2 2 2 2 0 1 1 0 1]
time: 101.36045880062589 s
16 bad initial partitions
>>> Reg
(array([0, 1, 2, 2, 1, 2, 2, 2, 2, 2, 0, 1, 1, 0, 1]),
11.0,
array([[0., 2., 8.],
       [0., 1., 0.],
       [0., 0., 0.]]),
array([[1, 1, 1],
       [1, 1, 3],
       [1, 3, 3]]))
>>> 
</code>

{{vlado:work:pics:classgreg11.png?350}}
{{vlado:work:pics:classmreg11.png?300}}

===== Structural BM =====

<code python>
>>> Types = np.array([{NUL,COM}]*(nC*nC),dtype='O').reshape(nC,nC)
>>> Rez = searchBM(500)
netsBM: 2020-08-15 00:12:05.699613
1 0 34.0 
  C : [0 2 1 2 2 0 0 2 0 1 1 2 0 1 0] 
  C*: [0 0 1 2 1 2 1 2 1 1 0 2 0 0 0]
6 2 28.0 
  C : [2 1 1 1 0 1 0 2 1 2 2 0 2 1 1] 
  C*: [1 1 0 2 1 2 0 1 0 0 1 2 1 1 1]
time: 127.35457458549575 s
2 bad initial partitions     
>>> Rez
(array([1, 1, 0, 2, 1, 2, 0, 1, 0, 0, 1, 2, 1, 1, 1]),
28.0,
array([[2., 6., 0.],
       [9., 4., 1.],
       [0., 6., 0.]]),
array([[2, 1, 1],
       [2, 1, 1],
       [1, 1, 2]]))
>>> 
</code>
{{vlado:work:pics:classgstr28.png?350}}
{{vlado:work:pics:classmstr28.png?300}}

====== To do ======

  - add the normal form of partition (done, August 15, 2020)
  - collect all different (locally) minimal solutions (done, August 15, 2020)
  - include undirected simple networks (done, August 15, 2020)


See [[.:gbm|Generalized blockmodeling based on Nets]].