非計量多次元尺度法
Nonmetric Multidimensional Scaling
非類似度データに対して、それに対応する対象の布置を与える座標
を求める。座標と距離の関係は、Minkowski distance (Borg & Groenen, 1997)、Minkowski-metric (Kruskal, 1964)と呼ばれている次式で与える。
いま、非類似度データに対応するとして与えられた距離を
、座標値から算出された距離を
で表すとき、これらの差の２乗和
は座標値の関数である。距離に対して関数（２）を最小にする座標値を極値探索法により求める。極値探索法として、scipy.optimize.minimizeのPowell法を用いた。
極値探索法は初期値が必要であるが、多次元尺度法の場合、初期値の影響を強く受けやすい。本ウェブサイトでは、Torgerson (1958)の尺度法（Torgerson scaling）と呼ばれている方法を用いて初期値を設定した（Borg, Groenen & Mair, 2013, p. 81）。すなわち、距離データに対して２重中心化（double centering: ）と固有分解（characteristic
decomposition）によって初期値を設定した（Borg & Groenen, 2005;
Carroll and Wish, 1974）。距離行列を
、座標値行列を
で表すとき、ユークリッド距離（2-metric）の下で導かれる次式で求められる（岡本安晴「心理学データ分析と測定」勁草書房、pp. 197-200）。
式（３）で与えられる布置を初期値とする極値探索法によって、式（２）の値を極小にする座標値を求めるが、１回目の探索でよい極小値が得られるとは限らない。得られた解にランダムな変動を加えて再度極値探索を行うということを繰り返して適切と思われる解を選ぶ。
さらに、非計量多次元尺度法では、非類似度データに対して設定される距離は、推定された座標値から計算される距離に対応して更新される。本ウェブサイトのプログラムでは、この更新をKruskal’s up-and-down-block algorithm
(Kruskal, 1964, pp. 127-128)と呼ばれている方法で行っている。このための関数をリストBAに示した。
得られた解とデータとの差の大きさは、ストレス公式１（stress formula 1: Borg & Groenen, 2005）と呼ばれている次式で算出している。
プログラムは、r-metricのrを指定して実行するものと、rが無限大であるsup metric用のものの２つを用意した。これらのプログラムのソースコードファイルとサンプルデータファイルは、圧縮ファイルfilesNM.zipとして用意した。Minkowski r-metricにおいて、rが有限値の場合（City metric、Euclid distanceなど）とｒが無限の場合（sup metric）とに分けて説明する。
Minkowski
r-metric (1 <= r < infinity)
プログラムのソースコードは、使用モジュールをリストR1に、本体はリストR2に示した。
リストR1のスクリプトには、１次元のベクトルの関数から、多次元の配列を構成して、式（２）の誤差を求める関数がクラス型linear_func内に定義されている。
def f(self, z):
x =
z.reshape((self.n_obj, self.n_dim))
t_dist = coord2dist( x, self.r )
err
= np.sum((self.dist - t_dist)**2)
return err
これにより、ベクトルzを多次元配列xに変換して極値探索法を用いている。
Minkowski r-metric（１）は、r=1のときはcity-block、Manhattan metricなどと呼ばれ、次式で表される。
図１.１で表される９個の点A～Iに対して、式（５）による距離を算出して、さらに単調変換を行い図1.２のデータを作成した。
図１.１
図1.２のデータは、図1.１の布置から作成された非類似度データを配列の下三角形の領域のみに設定したもので、この表データをCSVファイルとして保存したものを入力データとして読み込む。CSV形式のファイルは、＜岡本安晴「いまさら聞けないPythonでデータ分析」丸善出版、pp. 33-35, 101-102＞で説明している。
図1.２
リストR2のプログラムNonMetricMDS_r_metric.pyを実行すると、まず入力データファイル名の設定が求められる（図1.３）。
図1.３
図1.３では、図１.２のデータファイルDataCity.csvが設定されている。次に、出力用テキストファイル名の設定が求められるので、適当なテキストファイル名を設定する。図1.３では、Results.txtが設定されている。出力ファイル名の設定に続いて、r-metricのrの設定が求められる。City block metricでの多次元尺度法であるので、図1.３では、１が設定されている。このrの設定後、読み込んだデータの表示が行われ、続いて、多次元尺度法で用いる次元数の設定が求められる。図1.4では、図1.1の布置に対応して２が設定されている。生データのときは、いくつかの次元数で尺度構成を行ってみて、適切な解を選ぶことになる。
図1.4
次元数が設定されると計算が始まり、求められた布置による距離を縦軸に、入力データ値（で表されている）を横軸にとった散布図（Shepard diagram）が表示される（図1.5）。
図1.5
保存ボタンをクリックすると、表示グラフを保存することができる。閉じるボタンをクリックすると、Windowは閉じられる。
図1.5のフォームを閉じると、求められた解による点の布置が表示される（図1.6）。
図1.6
図1.6のフォームを閉じると、次に表示する布置の横軸の次元dim1と縦軸の次元dim2の設定が求められる。dim1に０以下の整数値を設定すると布置の表示終了となり、よりよい布置の可能性を探るために尺度構成を続けるかどうかを聞いてくる（図1.7）。
図1.7
図1.5のShepard diagramおよび図1.6に示されている布置が図1.１の布置と比べてより満足のいくものが得られることを期待して尺度構成を続けることにする。極値探索法の場合、探索の出発点（初期値）を変えると、異なる点（極小点）に収束する可能性がある。尺度構成を続けるため、Cを入力する。
よりよい布置を求めるために尺度構成を繰り返して、２回目（Solution 2）に得られた布置によるShepard diagramが図1.8である。
図1.8
このときの点の布置は図1.9のようになっている。座標軸の反転を除いてよく図1.1の布置に一致している。座標軸の反転は、式（５）のモデルによる距離が座標軸の反転に関して不変であることによる。
図1.9
プログラムを終了するときは、図1.7の状態においてQを入力すればよい。
Minkowski r-metricは、r=2のときはユークリッド距離と呼ばれ、次式で表される。
プログラムの使用法は、r=2と設定することを除いて、上のcity block metricの場合と同様である。
図1.1の布置に対してユークリッド距離を算出して単調変換を行い、図1.10のデータを得た。
図1.10
プログラムNonMetricMDS_r_metric.pyを起動して、図1.10のCSV形式のファイルDataEuclid.csvを入力ファイルとして指定する（図1.11）。
図1.11
出力用テキストファイル名を適当に設定して、r-metricのrを２とする。次元数は２を設定した。計算が終了すると、まず入力データ（で表されている）を横軸、求められた布置から算出された距離を縦軸とする散布図（Shepard diagram）が表示される（図1.12）。
図1.12
図1.12のフォームを閉じると、次に対象の布置が表示される（図1.13）。
図1.13
ユークリッド距離は、布置の回転および反転に関して不変であるので、図1.13の布置は図1.１の布置とよく対応していると言える。
Minkowski r-metricにおいて、
の場合は次式で与えられる。
距離（７）は、
-metric、sup metric (Kruskal, 1964)、dominance metric (Borg & Groenen, 1997)などと呼ばれている。式（７）の場合は、座標値の差が最大のもののみが距離に寄与するので、布置と距離の関係が弱い。このため、まずr=5の場合をsup metricの近似として用いて極値探索を行い、その後、sup metricによる極値探索を行っている（リスト２参照）。プログラムのソースコードはウェブサイトの後ろにリスト２として示した。極値探索がr=5のr-metricによるものとsup
metricによるものの２段構成になっていることを除けば、上のr-metricの場合（リスト１）と同じである。
図1.1の布置に対して、sup metricによる距離データを作成して単調変換を行い図2.1のデータを得た。
図2.1
プログラムNonMetricMDS_sup_metric.pyを起動して、入力データファイルとして図2.1のファイルDataSup.csvを設定し、出力用テキストファイル名を適当に設定する（図2.2）。
図2.2
次元数を設定すると（図2.2では２）、計算が始まり、極値探索が終了するとShepard diagramが表示される（図2.3）。
図2.3
図2.3のフォームを閉じると、次に布置が表示される（図2.4）。
図2.4
反転を除いて図1.1の布置によく対応しているようであるが、よりよい解の得られることを期待して尺度構成を繰り返す。
図2.5のShepard diagramを得る。
図2.5
このときの布置は、図2.6のようである。
図2.6
尺度構成の繰り返しはここまでとするときは、図1.7の画面でQを入力すればよい。
参考文献
Borg, L. & Groenen, P. (1997). Modern multidimensional scaling: Theory and applications. Springer.
Borg, I. & Groenen, P.J.F. (2005). Modern multidimensional scaling: Theory and application, 2nd ed. Springer.
Borg, I., Groenen, P. J. F. & Mair, P. (2013). Applied multidimensional scaling. New York: Springer.
Carroll, J, D. & Wish, M. (1974). Models and methods for three-way multidimensional scaling. In D. H. Krantz, A. C. Atkinson, R. D. Luce, & P. Suppes (Eds.). Contemporary developments in mathematical psychology, Vol. II.: Measurement, psychophysics, and neural information processing (pp. 57-105). Freeman.
Kruskal, J. B. (1964). Nonmetric multidimensional scaling: A numerical method. Psychometrika, 29, 115-129.
以下のファイルは圧縮ファイルfilesNM.zipにまとめられている。
リストBA Up-and-down-blocks algorithmのスクリプト例(BlockAverage_v1.py)
import copy
def nm_transf(dist):
"""
The up-and-down-blocks algorithm
I.
Borg & P. J. F. Groenen(2005). Modern multidimensional scaling, 2nd ed.,
(p. 206).
J.
B. Kruskal (1964). Nonmetric multidimensional scaling.
Psychometrika, 115-129 (pp. 126-128).
"""
a_dist = copy.copy(dist)
n = len(a_dist)
#
# Check
a_dist[0],...,a_dist[n-1] for monotonicity
#
Lpos = 0 # Start from a_dist[0]
Rpos = 0
transf_done = False
while transf_done == False:
#
# Check
a_dist[Lpos],...,a_dist[Rpos] for monotonicity
#
Lok
= False
Rok
= False
#
# Check the
right value a_dist[Rpos]
# If necessary,
adjust the values
#
if
Rpos >= n - 1:
Rok = True
elif
a_dist[Rpos] <= a_dist[Rpos+1]:
Rok = True
else:
v = a_dist[Rpos] * (Rpos - Lpos + 1) + a_dist[Rpos+1]
v /= Rpos - Lpos + 2
Rpos += 1
for i in range(Lpos, Rpos+1):
a_dist[i] = v
#
# Check
the left value a_dist[Lpos]
# If
necessary, adjust the values
#
if
Lpos <= 0:
Lok = True
elif
a_dist[Lpos-1] <= a_dist[Lpos]:
Lok = True
else:
v = a_dist[Lpos] * (Rpos - Lpos + 1) + a_dist[Lpos-1]
v /= (Rpos - Lpos + 2)
Lpos -= 1
for i in range(Lpos, Rpos+1):
a_dist[i] = v
if Lok and
Rok: # Monotonicity is satisfied. No adjustment
is done
Rpos += 1
Lpos = Rpos
if
Rpos > n-1:
transf_done = True
return a_dist
リストR1 リストR2で用いられるモジュール（mymodule.py）
import numpy as np
def coord2dist( x, r ):
"""
Calculate distances from coordinate x
by
r-metric
"""
n_obj = len(x)
n_dim = len(x[0])
dist = np.zeros((n_obj,
n_obj))
for i in range(n_obj):
for
j in range(i):
if j == i:
dist[i][i] = 0.0
else:
v = 0.0
for t in range(n_dim):
v
+= (np.abs(x[i][t] - x[j][t])) ** r
dist[i][j] = v ** (1.0 / r)
dist[j][i] = dist[i][j]
return dist
def coord2dist_sup( x ):
"""
Calculate distances from coordinate x
by sup-metric
"""
n_obj = len(x)
n_dim = len(x[0])
dist = np.zeros((n_obj,
n_obj))
for i in range(n_obj):
for
j in range(i):
if j == i:
dist[i][i] = 0.0
else:
v = 0.0
for t in range(n_dim):
if v < np.abs(x[i][t] - x[j][t]):
v = np.abs(x[i][t] - x[j][t])
dist[i][j] = v
dist[j][i] = dist[i][j]
return dist
class linear_func:
"""
Calculate sum of squared differences between
the target distances and derived r-metric
distances from flattened coordinate z
"""
def __init__(self, n_obj,
n_dim, dist, r):
self.n_obj = n_obj
self.n_dim = n_dim
self.dist = dist
self.r = r
def f(self, z):
x =
z.reshape((self.n_obj, self.n_dim))
t_dist = coord2dist( x, self.r )
err
= np.sum((self.dist - t_dist)**2)
return err
class linear_func_sup:
"""
Calculate sum of squared differences between
the target distances and derived sup-metric
distances from flattened coordinate z
"""
def __init__(self, n_obj,
n_dim, dist):
self.n_obj
= n_obj
self.n_dim = n_dim
self.dist = dist
def f(self, z):
x =
z.reshape((self.n_obj, self.n_dim))
t_dist = coord2dist_sup( x )
err
= np.sum((self.dist - t_dist)**2)
return err
リストR2 r-metricによるMDSスクリプト例（NonMetricMDS_r_metric.py）
import csv
import numpy as np
import scipy.optimize as so
import matplotlib.pyplot as plt
from mymodule import *
from BlockAverage_v1 import * # as ba
fin_name = input('Input data
file(*.csv) = ')
fout_name = input('Output file
name(*.txt) = ')
fout = open(fout_name, 'w')
fout.write('Input date
file...{}\n'.format(fin_name))
r = float(input('For r-metric, r = '))
fout.write('For r-metric, r =
{}\n'.format(r))
with open(fin_name, 'r') as fin:
raw_data = [v for v in
csv.reader(fin)]
obj_labels = raw_data[0][1:]
n_obj = len(raw_data) - 1
print('n_obj = ', n_obj)
dist = np.zeros((n_obj, n_obj))
for i in range(n_obj):
for j in range(i + 1):
if j
== i:
dist[i][i] = 0.0
else:
dist[i][j] = float(raw_data[i + 1][j + 1])
dist[j][i]
= dist[i][j]
print('dist...\n', dist)
fout.write('\nDistance Data...\n')
for i in range(n_obj):
fout.write('{0:
<5s}'.format(obj_labels[i]))
for j in range(i):
fout.write(' {0: >7.2f}'.format(dist[i][j]))
fout.write('\n')
n_dim = int(input('Number of
dimensions = '))
dist /= dist.max()
#
# Double
centering and characteristic decomposition
# for
initial configuration
#
NegD2 = (dist ** 2) * (-0.5)
J = np.zeros((n_obj, n_obj))
for i in range(n_obj):
J[i][i] = 1.0
J -= np.ones((n_obj, n_obj)) / n_obj
XXp = J @ NegD2 @ J
# Double centering
w, v = np.linalg.eigh(XXp) # Characteristic decomposition
idx = np.argsort(-w)
w = w[idx]
v = (v.T[idx]).T
X_init = np.zeros((n_obj, n_dim))
for t in range(n_dim):
if w[t] > 0.0:
for
i in range(n_obj):
X_init[i][t] = v[i][t] * (w[t] ** 0.5)
d_delta = np.empty((n_obj * (n_obj -
1)) // 2)
t_pos = 0
for i in range(1, n_obj):
for j in range(i):
d_delta[t_pos] = dist[i][j]
t_pos += 1
SltnID = 1
while True:
iter = 0
x1 = None
Stress_1 = None
dhat_t_1 = dist # Initial values
dhat_t_0 = None
while True:
print('\niter = ', iter)
func
= linear_func(n_obj, n_dim, dhat_t_1, r)
if
iter <= 0:
if SltnID <= 1:
x0 = X_init.reshape(n_obj * n_dim)
else:
x0 = X_init.reshape(n_obj * n_dim) + \
np.random.normal(scale = 0.1, size = n_obj * n_dim)
else:
x0 = x1.reshape(n_obj * n_dim) + \
np.random.normal(scale = 0.1, size = n_obj * n_dim)
iter
+= 1
res
= so.minimize(func.f, x0, method='Powell',
options={'maxiter':1000, 'disp':True})
print('res =\n', res)
x1 =
res.x.reshape((n_obj, n_dim))
means = x1.mean(axis = 0)
x1 =
x1 - means
#
#
Estimated distances by the optimized coordinates
#
d_hat = np.empty((n_obj * (n_obj - 1)) // 2)
i_pos = np.empty((n_obj * (n_obj - 1)) // 2)
j_pos = np.empty((n_obj * (n_obj - 1)) // 2)
t_pos = 0
for i in range(1,
n_obj):
for j in range(i):
v = 0.0
for t in range(n_dim):
v += (np.abs(x1[i][t] - x1[j][t])) ** r
d_hat[t_pos] = v ** (1.0 / r)
i_pos[t_pos] = i
j_pos[t_pos] = j
t_pos += 1
#
#
Transformation by up-and-down blocks algorithm
#
d_pair = np.array([(u, v, ipos, jpos) for u, v, ipos, jpos in
zip(d_delta, d_hat, i_pos, j_pos)],
dtype = ([('delta', float), ('dhat', float),
('ipos', int), ('jpos', int)]))
"""
dtype = ([('delta', np.float), ('dhat', np.float),
('ipos', np.int), ('jpos', np.int)]))
"""
idx
= np.argsort(d_pair, order = ['delta', 'dhat'])
dhat
= d_pair[idx]['dhat']
i_pair = d_pair[idx]['ipos']
j_pair = d_pair[idx]['jpos']
transf_d = nm_transf(dhat) # Transformation
dhat_t = transf_d
delta = d_pair[idx]['delta']
dhat_ave = np.zeros((n_obj, n_obj))
for
p in range((n_obj * (n_obj - 1)) // 2):
dhat_ave[i_pair[p]][j_pair[p]] = transf_d[p]
dhat_ave[j_pair[p]][i_pair[p]] = transf_d[p]
dhat_ave /= dhat_ave.max()
dhat_t_0 = dhat_t_1
# The previous values
dhat_t_1 = dhat_ave
# The new estiamted
distances
diff_d = np.sum((dhat_t_0 - dhat_t_1) ** 2)
diff_d /= (n_obj * (n_obj - 1)) /2
print('\ndiff_d = {}\n'.format(diff_d))
if
(diff_d < 1.0e-3):
print('\niter = {}\n'.format(iter))
v_num = np.sum((dhat - dhat_t) ** 2)
v_den = np.sum(dhat ** 2)
Stress_1 = (v_num / v_den) ** 0.5
break
print('\nSolution-{}'.format(SltnID))
print('\nStress-1 = ',
Stress_1)
print('\nX =')
fout.write('\nSolution-{}\n'.format(SltnID))
fout.write('\nStress-1 =
{}\n'.format(Stress_1))
fout.write('\nX =\n')
for i in range(n_obj):
print('{0: <5s}'.format(obj_labels[i]), end = '')
fout.write('{0: <5s}'.format(obj_labels[i]))
for
t in range(n_dim):
print(' {0: >7.3f}'.format(x1[i][t]), end = '')
fout.write(' {0: >7.3f}'.format(x1[i][t]))
print()
fout.write('\n')
#
# Shepard
diagram
#
plt.plot(d_delta, d_hat,
'bo')
plt.plot(delta, dhat_t, 'b')
plt.xlabel(r'$\delta_{i,j}$', fontsize = 14)
plt.ylabel(r'$\hat{d}_{i,j}$', fontsize = 14)
plt.title('Shepard diagram:
Solution-{}'.format(SltnID), fontsize = 16)
plt.tight_layout(h_pad = 5, w_pad
= 5)
plt.show()
#
#
Configuration of objects
#
dim1 = 0
dim2 = 1
if dim2 >= n_dim:
dim2
-= 1
while True:
plt.figure(figsize = (7,7,))
x_coord = x1.T[dim1]
y_coord = x1.T[dim2]
for
px, py, s in zip(x_coord, y_coord, obj_labels):
plt.text(px, py, s, fontsize = 14)
plt.plot(x_coord, y_coord, 'bo')
plt.xlabel('Dim-{}'.format(dim1+1), fontsize = 14)
plt.ylabel('Dim-{}'.format(dim2+1), fontsize = 14)
plt.title('Map: Solution-{}'.format(SltnID), fontsize = 16)
plt.tight_layout(h_pad = 5, w_pad = 5)
plt.show()
print('For dimension, put non-positive integer if you want to stop
drawing a map')
dim1
= int(input('dim1 = ')) - 1
if
dim1 < 0:
break
if
dim1 > n_dim - 1:
dim1 = n_dim - 1
dim2
= int(input('dim2 = ')) - 1
if
dim2 < 0:
dim2 = 0
elif
dim2 > n_dim - 1:
dim2 = n_dim - 1
ck = input('C(ontinue) or
Q(uit) ? ')
if ck[0] in {'Q', 'q'}:
break
SltnID += 1
fout.close()
print('\n{} was
saved.\n'.format(fout_name))
リストS1 リストS2のスクリプトで用いられるモジュール（mymodule_sup.py）
import numpy as np
def coord2dist( x, r ):
"""
Calculate distances from coordinate x
by
r-metric
"""
n_obj = len(x)
n_dim = len(x[0])
dist = np.zeros((n_obj,
n_obj))
for i in range(n_obj):
for
j in range(i):
if j == i:
dist[i][i] = 0.0
else:
v = 0.0
for t in range(n_dim):
v += (np.abs(x[i][t] - x[j][t])) ** r
dist[i][j] = v ** (1.0 / r)
dist[j][i] = dist[i][j]
return dist
def coord2dist_sup( x ):
"""
Calculate distances from coordinate x
by sup-metric
"""
n_obj = len(x)
n_dim = len(x[0])
dist = np.zeros((n_obj,
n_obj))
for i in range(n_obj):
for
j in range(i):
if j == i:
dist[i][i] = 0.0
else:
v = 0.0
for t in range(n_dim):
if v < np.abs(x[i][t] - x[j][t]):
v = np.abs(x[i][t] - x[j][t])
dist[i][j] = v
dist[j][i] = dist[i][j]
return dist
class linear_func:
"""
Calculate sum of squared differences between
the target distances and derived r-metric
distances from flattened coordinate z
"""
def __init__(self, n_obj,
n_dim, dist, r):
self.n_obj = n_obj
self.n_dim = n_dim
self.dist = dist
self.r = r
def f(self, z):
x =
z.reshape((self.n_obj, self.n_dim))
t_dist = coord2dist( x, self.r )
err
= np.sum((self.dist - t_dist)**2)
return err
class linear_func_sup:
"""
Calculate sum of squared differences between
the target distances and derived sup-metric
distances
from flattened coordinate z
"""
def __init__(self, n_obj,
n_dim, dist):
self.n_obj = n_obj
self.n_dim = n_dim
self.dist = dist
def f(self, z):
x =
z.reshape((self.n_obj, self.n_dim))
t_dist
= coord2dist_sup( x )
err
= np.sum((self.dist - t_dist)**2)
return err
リストS2 Sup距離のMDSスクリプト（NonMetricMDS_sup_metric.py）
import csv
import numpy as np
import scipy.optimize as so
import matplotlib.pyplot as plt
from mymodule_sup import *
from BlockAverage_v1 import *
fin_name = input('Input data
file(*.csv) = ')
fout_name = input('Output file
name(*.txt) = ')
fout = open(fout_name, 'w')
fout.write('Input date
file...{}\n'.format(fin_name))
r = 5
# r-metric is used to
estimate an approximate configuration
with open(fin_name, 'r') as fin:
raw_data = [v for v in
csv.reader(fin)]
obj_labels = raw_data[0][1:]
n_obj = len(raw_data) - 1
print('n_obj = ', n_obj)
dist = np.zeros((n_obj, n_obj))
for i in range(n_obj):
for j in range(i + 1):
if j
== i:
dist[i][i] = 0.0
else:
dist[i][j] = float(raw_data[i + 1][j + 1])
dist[j][i] = dist[i][j]
print('dist...\n', dist)
fout.write('\nDistance Data...\n')
for i in range(n_obj):
fout.write('{0:
<5s}'.format(obj_labels[i]))
for j in range(i):
fout.write(' {0: >7.2f}'.format(dist[i][j]))
fout.write('\n')
n_dim = int(input('Number of
dimensions = '))
dist /= dist.max()
#
# Double
centering and characteristic decomposition
# for
initial configuration
#
NegD2 = (dist ** 2) * (-0.5)
J = np.zeros((n_obj, n_obj))
for i in range(n_obj):
J[i][i] = 1.0
J -= np.ones((n_obj, n_obj)) / n_obj
XXp = J @ NegD2 @ J
# Double centering
w, v = np.linalg.eigh(XXp) # Characteristic decomposition
idx = np.argsort(-w)
w = w[idx]
v = (v.T[idx]).T
X_init = np.zeros((n_obj, n_dim))
for t in range(n_dim):
if w[t] > 0.0:
for
i in range(n_obj):
X_init[i][t] = v[i][t] * (w[t] ** 0.5)
x0 = np.zeros(n_obj * n_dim)
d_delta = np.empty((n_obj * (n_obj -
1)) // 2)
t_pos = 0
for i in range(1, n_obj):
for j in range(i):
d_delta[t_pos] = dist[i][j]
t_pos += 1
SltnID = 1
while True:
iter = 0
x1 = None
Stress_1 = None
dhat_t_1 = dist # Initial values
dhat_t_0 = None
while True:
print('\niter = ', iter)
func
= linear_func(n_obj, n_dim, dhat_t_1, r) # Funciton with r-metric
func_sup = linear_func_sup(n_obj, n_dim, dhat_t_1) # Function with sup-metric
if
iter <= 0:
if SltnID <= 1:
x0 = X_init.reshape(n_obj * n_dim)
else:
x0 = X_init.reshape(n_obj * n_dim) + \
np.random.normal(scale = 0.1, size = n_obj * n_dim)
else:
x0 = x1.reshape(n_obj * n_dim) + \
np.random.normal(scale = 0.1, size = n_obj * n_dim)
iter
+= 1
res
= so.minimize(func.f, x0, method='Powell',
options={'maxiter':1000, 'disp':True})
res
= so.minimize(func_sup.f, res.x, method='Powell',
options={'maxiter':1000, 'disp':True,
'xtol':1.0e-9, 'ftol':1.0e-14})
x1 =
res.x.reshape((n_obj, n_dim))
means = x1.mean(axis = 0)
x1 =
x1 - means
#
#
Estimated distances by the optimized coordinates
#
d_hat = np.empty((n_obj * (n_obj - 1)) // 2)
i_pos = np.empty((n_obj * (n_obj - 1)) //
2)
j_pos = np.empty((n_obj * (n_obj - 1)) // 2)
t_pos = 0
for
i in range(1, n_obj):
for j in range(i):
v = 0.0
for t in range(n_dim):
if
v < np.abs(x1[i][t] - x1[j][t]):
v = np.abs(x1[i][t] - x1[j][t])
d_hat[t_pos] = v
i_pos[t_pos] = i
j_pos[t_pos] = j
t_pos += 1
#
#
Transformation by up-and-down blocks algorithm
#
d_pair = np.array([(u, v, ipos, jpos) for u, v, ipos, jpos in
zip(d_delta, d_hat, i_pos, j_pos)],
dtype = ([('delta', float), ('dhat', float),
('ipos', int), ('jpos', int)]))
idx
= np.argsort(d_pair, order = ['delta', 'dhat'])
dhat
= d_pair[idx]['dhat']
i_pair = d_pair[idx]['ipos']
j_pair = d_pair[idx]['jpos']
transf_d = nm_transf(dhat)
dhat_t = transf_d
delta = d_pair[idx]['delta']
dhat_ave = np.zeros((n_obj, n_obj))
for
p in range((n_obj * (n_obj - 1)) // 2):
dhat_ave[i_pair[p]][j_pair[p]] = transf_d[p]
dhat_ave[j_pair[p]][i_pair[p]] = transf_d[p]
dhat_ave /= dhat_ave.max()
dhat_t_0 = dhat_t_1
dhat_t_1 = dhat_ave
#
The new estimated distances
diff_d = np.sum((dhat_t_0 - dhat_t_1) ** 2)
diff_d /= (n_obj * (n_obj - 1)) /2
if
(diff_d < 1.0e-3):
print('\niter = {}\n'.format(iter))
v_num = np.sum((dhat - dhat_t) ** 2)
v_den = np.sum(dhat ** 2)
Stress_1 = (v_num / v_den) ** 0.5
break
print('\nSolution-{}'.format(SltnID))
print('\nStress-1 = ',
Stress_1)
print('\nX =')
fout.write('\nSolution-{}\n'.format(SltnID))
fout.write('\nStress-1 =
{}\n'.format(Stress_1))
fout.write('\nX =\n')
for i in range(n_obj):
print('{0: <5s}'.format(obj_labels[i]), end = '')
fout.write('{0: <5s}'.format(obj_labels[i]))
for
t in range(n_dim):
print(' {0: >7.3f}'.format(x1[i][t]), end = '')
fout.write(' {0: >7.3f}'.format(x1[i][t]))
print()
fout.write('\n')
#
#
Shepard diagram
#
plt.plot(d_delta, d_hat,
'bo')
plt.plot(delta, dhat_t, 'b')
plt.xlabel(r'$\delta_{i,j}$', fontsize = 14)
plt.ylabel(r'$\hat{d}_{i,j}$', fontsize = 14)
plt.title('Shepard diagram:
Solution-{}'.format(SltnID), fontsize = 16)
plt.tight_layout(h_pad = 5,
w_pad = 5)
plt.show()
#
#
Configuration of objects
#
dim1 = 0
dim2 = 1
if dim2 >= n_dim:
dim2
-= 1
while True:
plt.figure(figsize = (7,7,))
x_coord = x1.T[dim1]
y_coord = x1.T[dim2]
for
px, py, s in zip(x_coord, y_coord, obj_labels):
plt.text(px, py, s, fontsize = 14)
plt.plot(x_coord, y_coord, 'bo')
plt.xlabel('Dim-{}'.format(dim1+1), fontsize = 14)
plt.ylabel('Dim-{}'.format(dim2+1), fontsize = 14)
plt.title('Map: Solution-{}'.format(SltnID), fontsize = 16)
plt.tight_layout(h_pad = 5, w_pad = 5)
plt.show()
print('For dimension, put non-positive integer if you want to stop
drawing a map')
dim1
= int(input('dim1 = ')) - 1
if
dim1 < 0:
break
if
dim1 > n_dim - 1:
dim1 = n_dim - 1
dim2
= int(input('dim2 = ')) - 1
if
dim2 < 0:
dim2 = 0
elif
dim2 > n_dim - 1:
dim2 = n_dim - 1
ck = input('C(ontinue) or
Q(uit) ? ')
if ck[0] in {'Q', 'q'}:
break
SltnID += 1
fout.close()
print('\n{} was
saved.\n'.format(fout_name))