Kamada, Kawai のグラフ描画アルゴリズム

はしがき

今回はフロイド法のところで言及したグラフ描画アルゴリズムについて書きます。
私はグラフ描画の世界に関しては何も知らない素人で、このアルゴリズムも古典的な有名どころらしきものを適当に選んだだけなのですが、このアルゴリズムが掲載されている An Algorithm for Drawing General Undirected Graphs は今日の時点で Google Scholar 調べで 1198 の記事から引用されているというとても凄い論文です。

アルゴリズム

アルゴリズムはとてもシンプルなもので、グラフの各頂点の現在の描画位置を引数にとるエネルギー関数を定義し、その極小点を Newton-Raphson 法で求めるだけです。
頂点を結ぶ辺は、頂点の位置が定まった後に直線によって描画されます。
よって（入力が平面グラフであっても）辺は交差している場合があります。

エネルギー関数の詳細と擬似コードは後で書く

C＋＋プログラムの一例です。
非常に長くなってしまいましたが、下手にソースを分割すると後で使いにくくなるおそれがあるので、過去に似たようなコードが出てきたような部分も含めて、全て一つのファイルに収めました。
また計算結果を出力をどうするか迷ったのですが、HTMLとJavaScriptのコードを出力させてブラウザに表示するようにしました。
一般的なブラウザさえあれば、余計なソフトを追加でインストールしなくても大丈夫です。

ちなみに、このアルゴリズムは定義より明らかに、非連結なグラフに対して適用すると挙動がおかしくなります。
しかしランダムグラフで作っている入力グラフは平気で非連結になるので、サンプルの出力がおかしい場合は何回かリロードしてみてください。

#include <cmath>
#include <cstdlib>
#include <ctime>
#include <iomanip>
#include <iostream>
using namespace std;

inline double roll() { return (double)rand() / RAND_MAX; }
const double MAX = pow(2.0, 64.0);
double** allocate_square_matrix(int num_element);

int main() {
	srand((unsigned int)time(0));
	
	const int num_vertices = 15;
	const int& N = num_vertices;
	
	double** adjacent = allocate_square_matrix(N);
	double**& A = adjacent;
	
	// make up an input graph (make a random graph)
	const double p = 0.2;
	const double max_length = 20;
	const double min_length = 1;
	for (int v = 0; v < N; v++) {
		A[v][v] = MAX;
		for (int u = v+1; u < N; u++) {
			double l = (roll() < p)
			? (max_length - min_length) * roll() + min_length : MAX;
			A[v][u] = l;
			A[u][v] = l;
		}
	}
	
	cerr << "The adjacent matrix of the input graph" << endl;
	cerr.setf(ios::fixed, ios::floatfield);
	cerr << setprecision(1);
	for (int v = 0; v < N; v++) {
		for (int u = 0; u < N; u++) {
			cerr << setw(4);
			if (A[v][u] == MAX)
				cerr << " MAX ";
			else
				cerr << A[v][u] << " ";
		}
		cerr << endl;
	}
	cerr << endl;

  // floyd's method
  double** distance = allocate_square_matrix(N);
  double**& D = distance;
  for (int v = 0; v < N; v++)
    for (int u = 0; u < N; u++)
      if (v == u)
        D[v][u] = 0.0;
      else
        D[v][u] = A[v][u];

  for (int k = 0; k < N; k++)
    for (int v = 0; v < N; v++)
      for (int u = 0; u < N; u++)
        if (D[v][k] + D[k][u] < D[v][u])
          D[v][u] = D[v][k] + D[k][u];

	cerr << "The distance matrix" << endl;
	cerr.setf(ios::fixed, ios::floatfield);
	cerr << setprecision(1);
	for (int v = 0; v < N; v++) {
		for (int u = 0; u < N; u++) {
			cerr << setw(4);
			if (D[v][u] == MAX)
				cerr << " MAX ";
			else
				cerr << D[v][u] << " ";
		}
		cerr << endl;
	}
	cerr << endl;

  // Kamada and Kawai's algorithm
  double** l = allocate_square_matrix(N);
  double** k = allocate_square_matrix(N);
  double L0 = 500;
  double K = 10;
  double max_distance = 0;
  for (int v = 0; v < N; v++)
    for (int u = 0; u < N; u++)
      if (D[v][u] > max_distance)
        max_distance = D[v][u];
  double L = L0 / max_distance;
  for (int v = 0; v < N; v++)
    for (int u = 0; u < N; u++) {
      l[v][u] = L * D[v][u];
      k[v][u] = K / D[v][u] / D[v][u];
    }

  double* x = new double[N];
  double* y = new double[N];
  for (int v = 0; v < N; v++) {
    x[v] = L0 * roll();
    y[v] = L0 * roll();
  }

	cerr.setf(ios::scientific, ios::floatfield);
	cerr << setprecision(1);
  const double epsilon = 0.001;
  const double& E = epsilon;
  while (1) {
    double max_delta = 0;
    double max_dEx = 0;
    double max_dEy = 0;
    int argmax_delta = -1;
    for (int v = 0; v < N; v++) {
      double dEx = 0;
      double dEy = 0;
      for (int u = 0; u < N; u++)
        if (u != v) {
          double d = (x[v] - x[u]) * (x[v] - x[u]) + 
            (y[v] - y[u]) * (y[v] - y[u]);
          d = sqrt(d);
          dEx += k[v][u] * ((x[v] - x[u]) - l[v][u] * (x[v] - x[u]) / d);
          dEy += k[v][u] * ((y[v] - y[u]) - l[v][u] * (y[v] - y[u]) / d);
        }
      double delta = sqrt(dEx * dEx + dEy * dEy);
      if (delta > max_delta) {
        max_delta = delta;
        max_dEx = dEx;
        max_dEy = dEy;
        argmax_delta = v;
      }
    }
    if (max_delta <= epsilon)
      break;
    //cerr << max_delta << endl;

    double delta = max_delta;
    double dEx = max_dEx;
    double dEy = max_dEy;
    int v = argmax_delta;
    while (delta > epsilon) {
      double dEx2 = 0;
      double dEy2 = 0;
      double dExy = 0;
      double dEyx = 0;
      for (int u = 0; u < N; u++)
        if (u != v) {
          double d = (x[v] - x[u]) * (x[v] - x[u]) + 
            (y[v] - y[u]) * (y[v] - y[u]);
          d *= sqrt(d);
          dEx2 += k[v][u] * (1 - l[v][u] * (y[v] - y[u]) * (y[v] - y[u]) / d);
          dExy += k[v][u] * l[v][u] * (y[v] - y[u]) * (x[v] - x[u]) / d;
          dEyx += k[v][u] * l[v][u] * (x[v] - x[u]) * (y[v] - y[u]) / d;
          dEy2 += k[v][u] * (1 - l[v][u] * (x[v] - x[u]) * (x[v] - x[u]) / d);
        }

      double a, b, c, d, det;
      det = dEx2 * dEy2 - dExy * dEyx;
      a = dEy2 / det;
      b = -dExy / det;
      c = -dEyx / det;
      d = dEx2 / det;
      x[v] += -dEx * a + -dEy * b;
      y[v] += -dEx * c + -dEy * d;

      dEx = 0;
      dEy = 0;
      for (int u = 0; u < N; u++)
        if (u != v) {
          double d = (x[v] - x[u]) * (x[v] - x[u]) + 
            (y[v] - y[u]) * (y[v] - y[u]);
          d = sqrt(d);
          dEx += k[v][u] * ((x[v] - x[u]) - l[v][u] * (x[v] - x[u]) / d);
          dEy += k[v][u] * ((y[v] - y[u]) - l[v][u] * (y[v] - y[u]) / d);
        }
      delta = sqrt(dEx * dEx + dEy * dEy);
      //cerr << "  " << delta << endl;
    }
  }

  // HACK
  for (int v = 0; v < N; v++) {
    x[v] += 50;
    y[v] += 50;
  }

	cout.setf(ios::fixed, ios::floatfield);
	cout << setprecision(3);
  cout << "<html>" << endl;
  cout << "  <head>" << endl;
  cout << "    <title>A graph</title>" << endl;
  cout << "    <script type='text/javascript'>" << endl;
  cout << "    onload = function() {" << endl;
  cout << "      var canvas = document.getElementById('graphcanvas');" << endl;
  cout << "      if (!canvas || !canvas.getContext)" << endl;
  cout << "        return false;" << endl;
  cout << "      var ctx = canvas.getContext('2d');" << endl;
  cout << "      ctx.fillStyle = '#ffffff';" << endl;
  cout << "      ctx.strokeStyle = '#000000';" << endl;
  cout << "      ctx.textAlign = 'center';" << endl;
  for (int v = 0; v < N; v++)
    for (int u = v+1; u < N; u++)
      if (A[v][u] < MAX) {
        cout << "      ctx.beginPath();" << endl;
        cout << "      ctx.moveTo(" << x[v] << ", " << y[v] << ");" << endl;
        cout << "      ctx.lineTo(" << x[u] << ", " << y[u] << ");" << endl;
        cout << "      ctx.closePath();" << endl;
        cout << "      ctx.stroke();" << endl;
      }
  for (int v = 0; v < N; v++) {
    cout << "      ctx.beginPath();" << endl;
    cout << "      ctx.arc(" << x[v] << ", " << y[v] 
         << ", 10, 0, Math.PI*2, true);" << endl;
    cout << "      ctx.fill();" << endl;
    cout << "      ctx.stroke();" << endl;
    cout << "      ctx.strokeText('" << v
         << "', " << x[v] << ", " << (y[v]+4) << ");" << endl; // HACK
  }
  cout << "    };" << endl;
  cout << "    </script>" << endl;
  cout << "  </head>" << endl;
  cout << "  <body>" << endl;
  cout << "    <canvas id='graphcanvas' width='600' height='600'></canvas>" << endl;
  cout << "  </body>" << endl;
  cout << "</html>" << endl;

  return 0;
}

double** allocate_square_matrix(int num_element) {
	double** matrix = new double*[num_element];
	matrix[0] = new double[num_element * num_element];
	for (int v = 1; v < num_element; v++)
		matrix[v] = matrix[v-1] + num_element;
  return matrix;
}