※時間がかかるので注意
In [8]:
%matplotlib inline
import numpy as np
from scipy.spatial.distance import euclidean as euc
import matplotlib.pyplot as plt
class Person:
def __init__(self, S, a, p=0.5):
self.S = S
self.a = a
self.p = p
def gather(self):
"""make person to participate the meeting.
"""
self.ideas = self.has_idea()
def has_idea(self):
"""a person has self.S ideas with self.a dimension.
"""
return list(np.random.rand(self.S, self.a))
def chose_idea(self, idea, idea2=None):
def nearness1(x, y, z):
"""calculate nearness of x for (y, z)
by calculating a linear combination.
"""
alpha = 1.
beta = 1.
return alpha*euc(x, y) + beta*euc(x, z)
def nearness2(x, y, z):
"""calculate nearness of x for (y, z)
by distance between x and the dividing point of (y, z) with t.
"""
# t > 0
# t <= 1: interior division
# t > 1: exterior division
t = 0.5
x, y, z = np.array(x), np.array(y), np.array(z)
return euc(t*(y-x) + (1.-t)*(z-x), (0., 0.))
if len(self.ideas) == 0:
return False
# return min(d) and its idea_id
if idea2 == None:
return min([(euc(vec, idea), idea_id) for idea_id, vec in enumerate(self.ideas)])
else:
return min([(nearness1(vec, idea, idea2), idea_id)
for idea_id, vec in enumerate(self.ideas)])
class Meeting:
"""Simulate a meeting with "simple3" situation.
Give keyword arguments:
K = 20 # Time limit
N = 6 # a number of participants
S = 10 # a number of ideas for each participants
a = 2 # the dimension of an idea
p = 0.5 # probability that a person speak
draw = True # draw image or don't
Output:
self.minutes: list of
( idea(which is vector with a dimension)
, who(person_id in the list "self.membes"))
self.k: stopped time (=len(self.minutes))
"""
def __init__(self, K=20, N=6, S=10, a=2, p=0.5, case=2):
self.K = K
self.N = N
self.S = S
self.a = a
self.p = p
self.case = case # case in the above cell: 2, 3, 4 or 5
if not self.case in [2, 3, 4, 5]:
raise ValueError
self.members = []
self.minutes = [] # list of (idea, who)
self.k = 0
def gather_people(self):
"""gather people for the meeting.
You can edit what ideas they have in here.
"""
for n in range(self.N):
person = Person(self.S, self.a, self.p)
# person.has_idea = some_function()
# some_function: return list of self.S arrays with dim self.a.
person.gather()
self.members.append(person)
self.members = np.array(self.members)
def progress(self):
"""meeting progress
"""
self.init()
preidea = self.subject
prepreidea = None
self.k = 1
while self.k < self.K + 1:
# l: (distance, speaker, idea_id) list for who can speak
l = []
for person_id, person in enumerate(self.members):
# chosed: (distance, idea_id)
chosed = person.chose_idea(preidea, prepreidea)
if chosed:
l.append((chosed[0], person_id, chosed[1]))
# if no one can speak: meeting ends.
if len(l) == 0:
print "no one can speak."
break
i = np.array([(person_id, idea_id)
for distance, person_id, idea_id in sorted(l)])
for person_id, idea_id in i:
rn = np.random.rand()
if rn < self.members[person_id].p:
idea = self.members[person_id].ideas.pop(idea_id)
self.minutes.append((idea, person_id))
if self.case == 3:
preidea = idea
elif self.case == 4:
prepreidea = idea
elif self.case == 5:
prepreidea = preidea
preidea = idea
self.k += 1
break
else:
self.minutes.append((self.subject, self.N))
self.k += 1
self.minutes = np.array(self.minutes)
def init(self):
self.gather_people()
self.subject = np.random.rand(self.a)
self.minutes.append((self.subject, self.N))
In [10]:
def myplot1(x, y, xfit=np.array([]), yfit=np.array([]), param=None,
scale=['linear', 'linear', 'log', 'log'], case=[2, 3, 4, 5]):
"""my plot function
x: {'label_x', xdata}
xdata: numpy array of array
y: {'label_y', ydata}
ydata: numpy array of array
param: {'a': 10, 'b': 20}
"""
if param:
s = [r'$%s = %f$' % (k, v) for k, v in param.items()]
label = s[0]
for _s in s[1:]:
label += ", " + _s
label_x, xdata = x.items()[0]
label_y, ydata = y.items()[0]
if len(scale)%2 == 1:
raise ValueError("'scale' must be even number")
fignum = len(scale)/2
figsize_y = 7 * fignum
fig = plt.figure(figsize=(10, figsize_y))
ax = []
for num in range(fignum):
ax.append(fig.add_subplot(fignum, 1, num+1))
for i, data in enumerate(zip(xdata, ydata)):
ax[num].plot(data[0], data[1], label="case: %d" % case[i])
if len(xfit):
ax[num].plot(xfit, yfit, label=label)
ax[num].legend(loc='best')
ax[num].set_xlabel(label_x)
ax[num].set_ylabel(label_y)
ax[num].set_xscale(scale[2*(num)])
ax[num].set_yscale(scale[2*(num)+1])
plt.show()
In [11]:
def calc_ave_dist_btw_nodes(case, trial, N):
_tmp = []
for t in range(trial):
meeting = Meeting(K=30, N=N, S=50, a=2, p=1.0, case=case)
meeting.progress()
tmp = 0
dist_btw_nodes= []
p0 = meeting.minutes[0]
for p1 in meeting.minutes[1:]:
if p1[1] != N:
dist_btw_nodes.append(euc(p0[0], p1[0]))
p0 = p1
_tmp.append(np.average(dist_btw_nodes))
_tmp = np.array(_tmp)
return np.average(_tmp)
def N_phi2(N, case=2, trial=100):
ave_dist_btw_nodes = [calc_ave_dist_btw_nodes(case, trial, _N) for _N in N]
return ave_dist_btw_nodes
In [17]:
import multiprocessing
case = [2, 3, 4, 5]
N = np.array([2**a for a in range(11)])
def wrapper(arg):
return arg[0](arg[1], arg[2], arg[3])
jobs = [(N_phi2, N, c, 100) for c in case]
process = multiprocessing.Pool(6)
ydata = np.array(process.map(wrapper, jobs))
In [18]:
myplot1({r'$N$': np.array([N]*len(case))}, {r'$\phi$':ydata}, case=case)
In [19]:
from scipy.optimize import leastsq
def myfit(fit_func, parameter, x, y, xmin, xmax, du=0, scale=['linear', 'linear', 'log', 'log']):
"""my fitting and plotting function.
fit_func: function (parameter(type:list), x)
parameter: list of tuples: [('param1', param1), ('param2', param2), ...]
x, y: dict
xmin, xmax: float
du: dataused --- index of the data used for fitting
"""
xkey, xdata = x.items()[0]
ykey, ydata = y.items()[0]
def fit(parameter, x, y):
return y - fit_func(parameter, x)
# use x : xmin < x < xmax
i = 0
while xdata[du][i] < xmin:
i += 1
imin, imax = i, i
while xdata[du][i] < xmax:
i += 1
imax = i
paramdata = [b for a, b in parameter]
paramkey = [a for a, b in parameter]
res = leastsq(fit, paramdata, args=(xdata[du][imin:imax], ydata[du][imin:imax]))
for i, p in enumerate(res[0]):
print parameter[i][0] + ": " + str(p)
fitted = fit_func(res[0], xdata[du][imin:imax])
fittedparam = dict([(k, v) for k, v in zip(paramkey, res[0])])
myplot1(x, y, xdata[du][imin:imax], fitted, param=fittedparam, scale=scale)
In [26]:
param = [('a', -0.5), ('b', 0.)]
xmin, xmax = 0., 1024
x = {r'$N$': np.array([N]*len(case))}
y = {r'$\phi$': ydata}
def fit_func(parameter, x):
a = parameter[0]
b = parameter[1]
return np.power(x, a)*np.power(10, b)
myfit(fit_func, param, x, y, xmin, xmax, du=2, scale=['log', 'log'])