2週間休み

明日からは平日の夕方以降がすべて某作業に忙殺されるため、学習は何も出来ないし休み時間もない。2週休みをうまく活用できたかと言われると疑問だが、まぁ休んだのなら休んだでいいのだろう。さて、この2週間は

・ポスター作った

・3連休だらだらした
・ポスター修正した

・UdemyのPyTorch講座をやった
・そこでjupyterの気持ちよさを知る
・ポスター修正した
・jupyterで何かしたい

・jupyterを使うために複雑ネットワーク解析関連を気分転換に実装しなおした。

PyTorch講座は、以前Chainer入門しようとしたが断念していたのでその埋め合わせ。結局のところ型どおりぶち込めばやってくれるっていうことですね、という程度の理解で視聴終了しました。
特にdeeplearningやる機会もなさそうなのでまた必要になったらやります。

で、jupyter欲を発散するためにネットワーク解析をまた実装した次第。
(後日追記)barabasi-albertモデルはググって出てきたものを実装したけど実際には違うらしいので修正します。

--------------------------------------------------------
#-*- coding:utf-8 -*-

import numpy as np
import matplotlib.pyplot as plt
import networkx as nx
import networkx as nx
import matplotlib.pyplot as plt
import pandas as pd
inf=np.float('inf')
global inf

def display_net(matrix):
    matrix=pathnet2connet(matrix)
    graph=nx.from_numpy_matrix(matrix)
    nx.draw(graph)
    plt.show()

def degree(matrix):
    matrix=pathnet2connet(matrix)
    return np.sum(matrix, 1)

def dykstra_matrix(matrix, start):
    matrix=connet2pathnet(matrix)
 
    end=np.zeros(matrix.shape[0])+inf
    pre_end=end.copy()
    end[start]=0
 
    while np.sum(end==pre_end)!=end.shape[0]:
        pre_end=end.copy()
        length=matrix+end.reshape(matrix.shape[0], 1)
        length[:, end!=inf]=inf
        if np.min(length)!=inf:
            end[np.argmin(length)%matrix.shape[0]]=np.min(length)
 
    return end

def arange(n):
    a=[]
    for i in range(n):
        a.append(i)
 
    return a

def list2matrix(list, n):
    list=np.array(list)
    matrix=np.zeros((n, n))
    for i in range(list.shape[0]):
        matrix[np.int(list[i, 0]), np.int(list[i, 1])]=1
 
    matrix=matrix+matrix.T
    return matrix

def matrix2list(matrix):
    matrix=pathnet2connet(matrix)
    list=[]
    for j in range(matrix.shape[0]):
        for i in range(j):
            if matrix[i, j]>0:
                list.append([i, j])
    return list

def pathnet2connet(matrix):
    matrix[matrix==inf]=0
    return matrix

def connet2pathnet(matrix):
    matrix[matrix==0]=inf
    for i in range(matrix.shape[0]):
        matrix[i, i]=0
 
    return matrix

def mk_WSmodel(n, k, p):
    #Watts-Strogatz modelのランダムネットワーク生成
    #input
    #n: ノード数
    #k: 次数
    #p: リンク張り替え確率
 
    matrix=mk_extended_cycle(n, k)
    matrix=change_link(matrix, p)
    matrix[matrix>0]=1
    return matrix

def mk_extended_cycle(n, k):
    #拡張サイクルモデルを作成する関数
    #input
    #n: ノード数
    #k: 次数(つまり前後にk/2本のリンクが張られる)
 
    #output: n*nの隣接行列
    matrix=np.zeros((n, n))
    k_half=np.int(k/2)
    for i in range(k_half):
        matrix[i, (n-k_half+i):]=1
        matrix[i, 0:i]=1
 
    i=k_half
 
    while i<n:
        matrix[i, (i-k_half):i]=1
        i=i+1
     
    matrix=matrix+matrix.T
    return matrix

def change_link(matrix, rate):
    n=matrix.shape[0]
    list=matrix2list(matrix)
    list=np.array(list)
 
    for i in range(list.shape[0]):
        edge=np.random.choice([0, 1], p=[0.5, 0.5])
        probability=[]
        for j in range(n):
            probability.append(rate/(n-2))
     
        probability[np.int(list[i, edge])]=1-rate
        probability[np.int(list[i, -(edge-1)])]=0
        num=arange(n)
        candidate=np.random.choice(num, p=probability)
        list[i, edge]=candidate
 
    matrix=list2matrix(list, n)
    return matrix

def cal_L(matrix):
    matrix=connet2pathnet(matrix)
    mean_length=np.zeros(matrix.shape[0])
    n=mean_length.shape[0]
 
    for i in range(n):
        mean_length[i]=np.sum(dykstra_matrix(matrix, i))/(n-1)
 
    L=np.mean(mean_length)
    return L

def cal_C(matrix):
    k=degree(matrix)
    t=[]
    for i in range(matrix.shape[0]):
        if k[i]>1:
            t.append(triangle(matrix, i))
        else:
            t.append(0)
            k[i]=2
 
    t=np.array(t)
    c=2*t/(k*(k-1))
    C=np.mean(c)
    return C

def triagle(matrix, i):
    counter=0
    matrix=pathnet2connet(matrix)
    n=matrix.shape[0]
    for j in range(n):
        for k in range(n):
            counter+=(matrix[i, j]*matrix[j, k]*matrix[k, i]>0)
    counter/=2
    return counter

def mk_random_network(n, l):
    list=[]
 
    i=2
    while i<n:
        counter=1
        while counter*i<=n:
            list.append(link_half2half((counter-1)*i, counter*i, i/2))
            counter=counter+1
     
        if (n-(counter-1)*i)>(i/2):
            list.append(link_half2half((counter-1)*i, n, i/2))
     
        i=i*2
 
    list.append(link_half2half(0, n, i/2))
 
    list=np.array(list)
    matrix=list2matrix(list, n)
 
    matrix=addlink(matrix, np.int(l-list.shape[0])   )   
     
    return matrix

def link_half2half(start, end, half):
    a=np.random.choice(np.arange(start, start+half))
    b=np.random.choice(np.arange(start+half, end))
    return [a, b]

def addlink(matrix, l):
    matrix=connet2pathnet(matrix)
    candidate=[]
    for j in range(matrix.shape[0]):
        for i in range(j):
            if matrix[i, j]==inf:
                candidate.append([i, j])
 
    candidate=np.array(candidate)
    list=np.random.choice(np.array(candidate).shape[0], l, replace=False)
    list=np.array(list)
 
    for i in range(list.shape[0]):
        x=candidate[np.int(list[i]), 0]
        y=candidate[np.int(list[i]), 1]
        matrix[x, y]=1
        matrix[y, x]=1
 
    return matrix

def smallworldness(matrix):
    count=100
    n=matrix.shape[0]
    matrix[matrix!=1]=0
    l=np.sum(matrix)/2
    S=(cal_C(matrix)/c_rand(n, l, count))/(cal_L(matrix)/l_rand(n, l, count))
    return S

def c_rand(n, l, count):
    c=0
    for i in range(count):
        c+=cal_C(mk_random_network(n, l))
 
    c/=count
    return c

def l_rand(n, l, count):
    l=0
    for i in range(count):
        l+=cal_L(mk_random_network(n, l))
 
    l/=count
    return l

def mk_BAmodel(m0, m, N):
    matrix=np.zeros((N, N))
    matrix[0:m0, 0:m0]=1
    for i in range(m0):
        matrix[i, i]=0

    i=m0:
    while i<N:
        probability=degree(matrix)
        probability=probability/np.sum(probability)
        list=np.random.choice(np.arange(N), m, p=probability)
        for j in range(m):
            matrix[i, list[j]]=1
            matrix[list[j], i]=1

        i=i+1

    return matirx

def mk_modifiedBAmodel(N):
    list=np.array([0, 1]).reshape(1, -1)
    n=2
 
    while n<N:
        matrix=list2matrix(list, n)
        d=degree(matrix)
        probability=d/np.sum(d)
        selection=np.random.choice(arange(n), p=probability)
        list=np.r_[list, np.array([n, selection]).reshape(1, -1)]
        n+=1
 
    matrix=list2matrix(list, n)
    return matrix
-----------------------------------------------
続いて検証。
Watts-Strogatzモデルのpを変化させた時のcharacteristic path lengthとclustering coefficientの推移をみる。
-----------------------------------------------
list=np.arange(0, 0.4, 0.05)
time=5
n=100
k=8
l=n*k/2
C_result=[]
L_result=[]

for p in list:
    C=0
    L=0
    for t in range(time):
        model=mk_WSmodel(n, k, p)
        C+=cal_C(model)/time
        L+=cal_L(model)/time
 
    C/=c_rand(n, l, time)
    L/=l_rand(n, l, time)
    C_result.append(C)
    L_result.append(L)

df1=pd.DataFrame({'p-value': list, 'clustering coefficient': C_result})
df2=pd.DataFrame({'p-value': list, 'characteristic path length': L_result})
plt.plot('p-value', 'clustering coefficient', data=df1)
plt.plot('p-value', 'characteristic path length', data=df2)
plt.show()

-----------------------------------------------

合ってそう。
次修正Barabasi-Albertモデル。アルゴリズム的にclustering coefficientはゼロなので、まぁグラフ見てなんとなく合っていればとりあえずいいか。
-----------------------------------------------
display_net(mk_modifiedBAmodel(40))
-----------------------------------------------

良さげ。
次数分布と次数相関係数も求めてみる。
-----------------------------------------------
def degree_distribution(matrix):
    deg=degree(matrix)
    result=np.zeros(np.int(np.max(deg)))
    for i in range(result.shape[0]):
        result[i]=np.sum(deg==(i+1))
 
    return result

def assortativity(matrix):
    #assortativity公式を(A-C)/(B-C)と見立ててそれぞれ求める
    #前準備
    matrix=pathnet2connet(matrix)
    N=matrix.shape[0]
    M=np.sum(matrix)
    A=0
    B=0
    C=0
 
    #先にdegreeを求めておく
    degree_vector=degree(matrix)
 
    #まずAから
    A=np.sum(np.dot(degree_vector, degree_vector.T)*matrix)/(4*M)
 
    #次B
    #i, jノードのそれぞれの2乗の和を(i,j)要素にもつ行列sum_sqを作成
    sq_degree=degree_vector*degree_vector
    sum_sq=ij_sum(sq_degree)
    B=np.sum(sum_sq*matrix)/(2*M)
 
    #次C
    ijsum=ij_sum(degree_vector)
    C=(np.sum(ijsum*matrix)/(2*M))**2
 
    assortativity=(A-C)/(B-C)
    return assortativity
 
#引数ベクトルのi番目、j番目の要素の和を(i, j)に持つ行列を返す
def ij_sum(vector):
    n=vector.shape[0]
    ij_sum=np.zeros((n, n))
 
    for i in range(n):
        for j in range(n):
            ij_sum[i, j]=vector[i]+vector[j]
 
    return ij_sum
-----------------------------------------------
assortatibityに関して、リンク数は重複カウントするのかどうかが微妙。とりあえず全部カウントすることにして実装した。真偽は定かでない。
-----------------------------------------------
N=300
total=np.zeros(N-1)
for i in range(100):
    model=mk_modifiedBAmodel(N)
    result=degree_distribution(model)
    for j in range(result.shape[0]):
        total[j]+=result[j]
 

total/=100
list=np.arange(299)


df=pd.DataFrame({'degree': list[0:10], 'count': total[0:10]})
plt.plot('degree', 'count', data=df)
plt.show()
assortativity(model)
-----------------------------------------------

5.2548947026591675