D言語 Advent Calendar 2012、23日目の記事でございます。今日はクリスマスイブの前日。なので、それにちなんで、巡回サラリーマン問題ならぬ「巡回サンタクロース問題」をD言語＋遺伝的アルゴリズムで解こうとおもいます。
サンタクロースの元ネタになったのは、カトリックの聖人ニコラオス。貧しい家の子どもを救うため、或いは冤罪被害者を救うため、彼にはクリスマスだけと言わずに、年がら年中奔走して貰いたいものです……しなさい。して下さい、ｵﾅｼｬｽ。
ということで、世界都市247の首都を最短経路＋αの条件で巡るプログラムを、D言語で作ってみたいと思います。
条件は、以下のとおりとします。
1 サンタの出発地、及び、最終到着地はグリーンランドとし、全世界の首都を通過する。
2 データはOpenGeoCodeから拝借した経度緯度情報を使用する。
※南極大陸など、人の居ない（場合によっては国家ですら無い）場所も通過する。
※中国、台湾、香港、マカオはそれぞれ別の国と看做される。ここに当研の政治的意図はない。
※イスラム圏国家も通過する。聖人が差別しちゃいけない。また、彼に地対空攻撃を行う不届き者もいない。
3 サンタの飛行高度は海面高度10000m、移動速度は毎年のようにサンタを追跡することに定評のあるNORADの公式発表である「新幹線の100倍の速度」=マッハ約23とする。
4 サンタは地球の自転に引きずられない。つまり、その場に浮いていると、徐々に西に流されていき、24時間後に同じ位置に戻ってくる。
5 空気抵抗、特にジェット気流は考慮しない。聖人は絶対空気抵抗なんかに負けたりしない!!
実装は、以下のとおりとします。
1 使用する探索法は、遺伝的アルゴリズムを使用した近似値探索とする（よって、最適解を期待しない）。
1-1 1個体の遺伝子数は247とし、それぞれに0～65534の数値と国コードをする。この数値が小さい順に通過する。特別措置として、グリーンランド（GRL）だけは65535を振り、必ず最後に通過するようにする。
1-2 適用度は「移動に要する時間の総計（秒）の負数」とし、高い（＝所要時間が小さい）サンタが生き残るようにする。
1-3 細かいパラメータは以下のとおり。
総計算回数 = 10000回
個体数 = 384体
両親選択 = 適用度トップ32体をルーレット方式で選択
交叉方法 = 1遺伝子ごとに、どちらを取得するかをランダムで選択
突然変異率 = 交叉時に、1遺伝子につき0.5%の確率で発生
ソースコード全文はhttp://dl.dropbox.com/u/1128855/tsp.zipから持って行って下さい。
では、ソースコードを見ていきましょう。
module cow;
import
ixa.geo.coord ,
ixa.geo.radian ,
ixa.geo.vector ,
std.algorithm ,
std.array ,
std.math ,
std.typecons ;
/// 地球計算条件
enum Earth : real
{
equator = 40000 * 1000 , // 赤道長さ[m]
radius = equator / 2 / PI , // 赤道半径[m]
jitensokudo = 2.0 * PI / (24*60*60) // 自転の角速度[π/s]
}
/// サンタ計算条件
enum SantaState : real
{
altitude = 10000 , // 飛行高度[m];
speed = 250.0 * 100 / 3.6 , // 移動速度( 新幹線の速度[km/h] * 100 のm/s換算)
sprad = (2*PI*speed) / (2*PI*(Earth.radius + altitude)) , // 角速度[π/s]
}
/// 国別コード
enum CountryCode
{
unknown = -1,
ABW = 0 ,
// 中略＼(^o^)／
ZWE = 246 ,
}
/// 国別コードの配列化
enum CountryCodes = {
CountryCode[] res;
foreach( name ; __traits(allMembers,CountryCode) ){
if( name.length == 3 ){
res ~= __traits(getMember,CountryCode,name);
}
}
return res;
}();
/// 国数
enum numCountry = { return CountryCodes.reduce!"a>=b?a:b" + 1; }();
struct CountryData
{
CountryCode _code; /// 国コード
string _name; /// 名称
string _capitalName; /// 首都名
uint _population; /// 人口
float _latitude; /// 緯度( +:北緯 -:南緯 )［60分法］
float _longtitude; /// 経度( +:東経 -:西経 )［60分法］
/// 地球核からその都市への正規化されたベクトルの算出
Vector!(float,3) toSurface()( auto ref Radian!float latitudeBias // 緯度バイアス値
) const
{
// ↓ベクトル軸(0,1,0)を時計回りにlatだけ回転させる変換行列を生成する
auto lat = degree(_latitude) + latitudeBias;
auto rLat = vector!float( 0,1,0).rotateMatrix( lat );
// ↓ベクトル軸(-1,0,0)を時計回りにlatだけ回転させる変換行列を生成する
auto lon = degree(_longtitude);
auto rLon = vector!float(-1,0,0).rotateMatrix( lon );
return rLat * rLon * coord!float(0,0,1) - Coord!(float,3).O ;
}
}
/// 国別データ
enum countries =
{
CountryData[CountryCode] r;
auto cs = [
CountryData( CountryCode.ABW,"Aruba","Oranjestad",107635,12.5,-69.96666667 ),
// 中略＼(^o^)／
CountryData( CountryCode.ZWE,"Zimbabwe","Harare",12619600,-20,30 )
];
foreach( c ; cs ){
r[c._code] = c;
}
return r;
}();
unittest {
static assert( countries.length == numCountry );
}
/// 移動元の国と移動先の国の開き角度を余弦で取得する（近いなら1、遠いなら-1）。
float distNsa( CountryCode to, // 移動先
CountryCode from, // 移動元
float nsec // 経過秒数
)
{
return countries[to].toSurface( Radian!float(nsec * Earth.jitensokudo) ) *
countries[from].toSurface( Radian!float(0) );
}
unittest {
// 30秒後のトルコと日本の距離の違いを算出
assert( distNsa( CountryCode.JPN, CountryCode.TUR, 30 ) < // = -0.224185 トルコ→日本は、自転の影響で遠くなる
distNsa( CountryCode.TUR, CountryCode.JPN, 30 ) // = -0.224046 日本→トルコは、逆に近くなる
);
}
/// 国fromから国toへの移動所要時間を計算する。計算結果は静的キャッシュされる。
float requiredTime( CountryCode to, CountryCode from )
{
static float[numCountry][numCountry] store;
if( store[to][from].isNaN ) {
auto A = (float N) => distNsa( to, from, N );
auto B = (float N) => cos(SantaState.sprad * N);
float refCalc( float ns, float r, int i ) {
auto n = pow( 2.0f, ns );
if( i > 0 ){
auto a = A(n);
auto b = B(n);
if( a < b ){
// 距離足らず
return refCalc( ns + r, r / 2, i - 1 );
}
else if( a > b ){
// 距離過多
return refCalc( ns - r, r / 2, i - 1 );
}
else {
return n;
}
} else {
return n;
}
}
store[to][from] = refCalc( 0, 10, 15 );
}
return store[to][from];
}
国別データの定義と、それにまつわる計算処理を行なっています。enum定義の際に、
enum SYMBOL = { return VALUE; }();
↑のような書き方をしていますが、CTFE（Compile Time Function Execution）によってコンパイル中に{}内の処理が行われ、VALUEが静的な値としてSYMBOLに格納されます。
次に、mainproc.dを見てみましょう。
/*(中略)*/
/// 個体
struct Santa
{
alias typeof(this) This;
ushort[247] _data; /// 遺伝子データ
/*(中略)*/
/// 遺伝情報から経路順序を動的配列で生成
CountryCode[] calcPassWay()
{
return _data[].zip( CountryCodes )
.array
.sort!((ref a,ref b) => a[0] < b[0])
.map!( a => a[1] )
.array;
}
/*(中略)*/
/// 交叉による自個体の更新
void crossover( ref This as, ref This bs )
{
struct S {
ushort[] a,b,c;
int opApply( scope int delegate(ref ushort,ref ushort,ref ushort) dg )
{
int r;
foreach( i ; 0 .. a.length ){
r = dg( a[i], b[i], c[i] );
if( r ) return r;
}
return r;
}
}
foreach( ref r,a,b ; S( _data, as._data, bs._data) ) {
if( uniform(0,2) ){
r = a;
} else {
r = b;
}
if( uniform(0,1000) < 5 ){
// 突然変異処理
r = cast(ushort)uniform( 0, 65535 );
}
}
_data[CountryCode.GRL] = 65535; // グリーンランドは必ず最後にする
}
}
ここでは、main関数とGAで処理する個体情報構造体Santaを定義しています。
遺伝情報から経路を生成するcalcPassWayでは、「遺伝情報１つことに国コードをつけたtupleにする(zip)」「遺伝情報の昇順に並べる(sort)」「国コードだけにする」という処理をしています。 UFCS導入前ですと、この処理は
return array(map!( a => a[1] )(sort!((ref a,ref b) => a[0] < b[0])(array(zip( _data[], CountryCodes )))));
等と書かねばならず、「処理の順序と記述の順序が逆」「カッコが多くて見づらい」「処理順に記述したかったら変数を用意しないといけない」など、マンモス哀れな記述を強いられていたんだ！
しかし、UFCSの導入によって、これらの記述的問題が解消されています。読みやすいコードで実装できるというのは、それだけでバグが減るものなのです。
最後の1つ、ga.dは割愛します。何のへんてつもないGAフレームワークなので……。
ということで、実行してみましょう。該当フォルダから、以下のコマンドをコマンドライン上から入れましょう。
rdmd mainproc
dmdの付属アプリrdmdは、ソースコードの依存関係も追跡してソースをコンパイル、実行します。とりあえずunittestを走らせたい、という時に重宝します。
結果は、以下の様な感じになります。
calculating 0 - 49 / 10000 : max = nan
calculating 50 - 99 / 10000 : max = -405455(中略)
calculating 9900 - 9949 / 10000 : max = -45162.8
calculating 9950 - 9999 / 10000 : max = -45162.8
-45153.4
CountryData(SPM, "Saint Pierre and Miquelon", "Saint-Pierre", 5831, 46.8333, -56.3333)
CountryData(BMU, "Bermuda", "Hamilton", 69080, 32.3333, -64.75)
CountryData(AIA, "Anguilla", "The Valley", 15423, 18.25, -63.1667)
CountryData(SXM, "Sint Maarten (Dutch part)", "Philipsburg", 37429, 18.0417, -63.0667)
CountryData(ATG, "Antigua and Barbuda", "Saint John's", 89018, 17.05, -61.8)(中略)
CountryData(AND, "Andorra", "Andorra la Vella", 85082, 42.5, 1.5)
CountryData(FRA, "France", "Paris", 65630692, 46, 2)
CountryData(JEY, "Jersey", "Saint Helier", 94949, 49.25, -2.16667)
CountryData(IRL, "Ireland", "Dublin", 4722028, 53, -8)
CountryData(GRL, "Greenland", "Nuuk", 57695, 72, -40)
最終適応度は-45153.4、即ち全世界を巡るのに45153秒、約12.5時間かかる、という結果が出ました。頑張ってください、サンタさん。
また、遺伝的アルゴリズムの悲しい習性として、結果は一定じゃありません。最悪だと、40%ぐらいの違いが出るのにはご注意下さい。
前の方々の記事を見て頂ければお分かりになりますが、D言語はまだ多くの問題を抱えています。それでも尚、不肖がD言語を選択する理由はやはり、「極めて高い生産性」です。twitterで不肖をフォローしている方はご存知かもしれませんが、2日目前までは、不肖はこの経路探索をA*アルゴリズムでやろうとしていました……が、収束までに予想以上の時間がかかるため、急遽遺伝的アルゴリズムでの探索に切り替えました。設計？哲学性？なにそれおいしいの？状態のカウボーイスタイルで大慌てで組んだのですが、結果はご覧のとおりです……できちゃいました。勿論、ちゃんとした設計手順を踏めば、更なる効率が望めるでしょう。そんなしょっぱい裏事情にまで対応しきる所に、D言語の良さがあるのだと不肖は思います。
本当なら、SDL+OpenGLを使った結果のビジュアライゼーション化、マルチスレッドによる高速化などもやってみたかったのですが、そこは不肖の未熟たるが故、またの機会とさせていただきたくご容赦願います。
ちなみに、ソース中に登場する「ixa」というモジュールは、不肖謹製のソースコード資源「Inexir eXtensional source Asset」です。無茶な使い方をするとすぐヘソを曲げるので、参考程度に御覧ください。Boost Licenseです。近々、githubにも展開予定です。
ということで、次は@psihyさんです。