【忘却のJava#4】組み合わせ最適化問題を解くプログラムをJavaで書く（１）

　今回は大学時代に研究書いていたC言語のコードをJavaで書いてみたいと思います。
その前にC言語で当時のプログラムを動かすところまで事前準備でやっておきます。
実際にJavaコード書くのは次回の予定です。

組み合わせ最適化問題とは？

ある問題において、ある制約の下で選べる組み合わせの中で、最良の組み合わせを選ぶ問題です。
例えば、大きさや重さといった制約のもとナップサックに荷物を詰めるナップサック問題があります。

【参考】
「組合せ最適化問題とは」

最小全域木問題 (Minimum Spanning Tree Problem)

組み合わせ最適化問題の一つとして最小全域木問題があります。
全域木とは与えられた重み付き無向グラフ中の頂点全てを含むような木のことです。
その全域木の中でも各辺の重みの和が最小となるものが最小全域木です。

【参考】
全域木とは：全域木 - Wikipedia
最小全域木とは：例えば数のようなものが最小全域木です。

f:id:wingfair:20190410224004p:plain — 左：重み付き無向グラフ　右：最小全域木

Kruskalのアルゴリズム

最小全域木の最適解を導くことができるアルゴリズムのうちの1つです。
その時点で最も重みの小さい辺を閉路ができないように追加していき、部分木を順次拡大させていくことで最終的には最小全域木になるといったものです。
詳しくは以下の【参考】をみてみてください。

【参考】
クラスカル法 - Wikipedia

これを題材にしてC言語で書いたものをJavaで書いてみようと思ってます。
今回はその前に昔を思い出しながらCで書いたプログラムを実行するところまでやってみます。

事前準備

C言語のプログラムのコンパイラであるgccをインストールしておきます。

# yum install gcc

Cのコード

久しぶりにCのコード見ましたが、printf()とかputs()が懐かしい。。
あとQuickSortをちゃんと書いてるのとか軽く感動。。

Kruskalの中ではUnion-Findというロジックを使ってます。
面白い手法なので興味ある人は以下の方のサイトとか参考にしてみてください。
【参考】
最小全域木（クラスカル法とUnionFind） - アルゴリズム講習会

graph.c

#include<stdio.h>
#include<stdlib.h>
#include"graph.h"

#define min(x,y)  ( x < y ? x : y )
#define max(x,y)  ( x > y ? x : y )

/* rerutn empty graph(node number = n) */
Graph EmptyGraph(int n)
{
  Graph g;
  int i, j;
  g.ord = n;
  g.adj = (int**)malloc(sizeof(int*) * n);
  for (i = 0; i < n; ++i) {
      g.adj[i] = (int*)malloc(sizeof(int) * n);
  }
  for (i = 0; i < n; i++) {
    for (j = 0; j < n; j++) {
      g.adj[i][j] = 0;
      g.adj[j][i] = 0;
    }
  }
  return g;
}

/* free graph array */
void FreeGraphArray(Graph g)
{
  int i;
  for (i = 0; i < g.ord; i++) {
    free(g.adj[i]);
    g.adj[i] = NULL;
  }
  free(g.adj);
  g.adj = NULL;
}

/* make graph from m[]  */
Graph ArraytoGraph(Graph g, int (*m)[g.ord])
{
  int i, j;
  for (i = 0; i < g.ord; i++) {
    for (j = 0; j < g.ord; j++) {
      g.adj[i][j] = m[i][j];
    }
  }
  return g;
}

/* Print Graph g */
void PrintWeightedGraph(Graph g)
{
  int i, j;
  for (i = 0; i < g.ord; i++) {
    for (j = 0; j < g.ord; j++) {
      printf("%4d", g.adj[i][j]);
    }
    puts("");
  }
}

/* return total edge number of Graph g */
int CountEdge(Graph g)
{
  int i, j;
  int s = 0;
  for (i = 0; i < g.ord; i++) {
    for (j = i + 1; j < g.ord; j++) {
      if (g.adj[i][j] > 0) {
        s++;
      }
    }
  }
  return s;
}

/* return parent of x */
Group *FindParent(Group *x)
{
  if (x->parent != x) {
    x->parent = FindParent(x->parent);
  }
  return x->parent;
}

/* marge x and y */
void Union(Group *x, Group *y)
{
  if (x->rank > y->rank) {
    y->parent = x;
    y->rank = x->rank;
  } else {
    x->parent = y;
    x->rank = y->rank;
    if (x->rank == y->rank) {
      y->rank++;
      x->rank++;
    }
  }
}

/* swap x and y */
void swap(Edge *x, Edge *y)
{
  Edge temp;
  temp = *x;
  *x = *y;
  *y = temp;
}

/* Quick Sort */
void QuickSortEdge(Edge *e, int low, int high)
{
  int i, j, mid;
  double pivot;

  mid = (low + high) / 2;
  i = low;
  j = high;
  pivot = e[mid].w;
  while(1) {
    while(e[i].w < pivot) { i++; }
    while(e[j].w > pivot) { j--; }
    if (i >= j) { break; }
    swap(&e[i++], &e[j--]);
  }
  if (low < i-1) { QuickSortEdge(e, low, i-1); }
  if (j+1 < high) { QuickSortEdge(e, j+1, high); }
}

/* make minimum spanning tree */
Graph Kruskal(Graph g, Graph tree)
{
  int i, j, k, a, b, edge_num, cn_comp = 0;
  Group *flm;
  Edge *ed;

  edge_num = CountEdge(g);
  //Group *flagment;
  ed = (Edge *)calloc(edge_num, sizeof(Edge));
  // add edge in set
  k = 0;
  for (i = 0; i < g.ord; i++) {
    for (j = i+1; j < g.ord; j++) {
      if (g.adj[i][j] > 0) {
        ed[k].w = MakeWeightValue(g.adj[i][j], i, j);
        ed[k].p1 = i;
        ed[k].p2 = j;
        k++;
      }
    }
  }
  flm = (Group *)malloc(sizeof(Group) * g.ord);
  for (i = 0; i < g.ord; i++) {
    flm[i].parent = &flm[i];
    flm[i].rank = 0;
    flm[i].id = i;
  }
  QuickSortEdge(ed, 0, edge_num-1);
  k = 0;
  while (cn_comp < g.ord - 1) {
    a = ed[k].p1;
    b = ed[k].p2;
    // find if cycle exists or not
    if (FindParent(&flm[a]) != FindParent(&flm[b])) {
      // add edge in MST
      tree.adj[a][b] = g.adj[a][b];
      tree.adj[b][a] = g.adj[b][a];

      // marge set
      Union(flm[a].parent, flm[b].parent);
      cn_comp++;
    }
    k++;
  }
  free(flm);
  flm = NULL;
  free(ed);
  ed = NULL;
  return tree;
}

/* return sum of tree's edge consts */
int TreeTotalEdgeCost(Graph g)
{
  int i, j, cost, sum = 0;
  for (i = 0; i < g.ord; i++) {
    for (j = i + 1; j < g.ord; j++) {
      cost = g.adj[i][j];
      if (cost > 0) {
        sum += cost;
      }
    }
  }
  return sum;
}

/* return no duplicate weight value */
int MakeWeightValue(int w, int id, int p)
{
  int value;
  int a, b, c;
  a = w * 100000;
  b = max(id, p) * 1000;
  c = min(id, p);
  value = a + b + c;

  return(value);
}

main_mst.c

#define N 6 // graph node number

#include<stdio.h>
#include "graph.h"

int main(void)
{
  Graph g, tree; // graph
  int tc; // total cost

  // make graph
  g = EmptyGraph(N);
  tree = EmptyGraph(N);

  // test data
  int maze[N][N] = {
    { 0, 7,15, 0, 0, 0 },
    { 7, 0,12, 6, 0, 0 },
    {15,12, 0,18,17, 2 },
    { 0, 6,18, 0,13, 0 },
    { 0, 0,17,13, 0,24 },
    { 0, 0, 2, 0,24, 0 },
   };

  // show test data
  g = ArraytoGraph(g, maze);
  puts("Graph (problem)");
  PrintWeightedGraph(g);

  // make MST by Kruskal Method
  tree = Kruskal(g, tree);
  puts("");
  puts("Tree (result)");
  PrintWeightedGraph(tree);

  // show result
  tc = TreeTotalEdgeCost(tree);
  printf("total cost:%d\n", tc);

  // close graph
  FreeGraphArray(g);
  FreeGraphArray(tree);

  return(0);
}

make

# make
gcc -g -Wall -O2 -c main_mst.c
gcc -o mst main_mst.o graph.o

実行結果

こんな感じで最適解を求めることができます。
この"Graph"問題の最小全域木は"Tree"で最小の重みは40という結果です。

# ./mst
Graph (problem)
   0   7  15   0   0   0
   7   0  12   6   0   0
  15  12   0  18  17   2
   0   6  18   0  13   0
   0   0  17  13   0  24
   0   0   2   0  24   0

Tree (result)
   0   7   0   0   0   0
   7   0  12   6   0   0
   0  12   0   0   0   2
   0   6   0   0  13   0
   0   0   0  13   0   0
   0   0   2   0   0   0
total cost:40

以上。

ailes blog