　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から持って行って下さい。
　では、ソースコードを見ていきましょう。

cow.d

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を見てみましょう。

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さんです。