Help us understand the problem. What is going on with this article?

D言語実戦編・遺伝的アルゴリズムで巡回サンタ問題を解く

More than 5 years have passed since last update.

 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では、「遺伝情報1つことに国コードをつけた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さんです。

Why not register and get more from Qiita?
  1. We will deliver articles that match you
    By following users and tags, you can catch up information on technical fields that you are interested in as a whole
  2. you can read useful information later efficiently
    By "stocking" the articles you like, you can search right away
Comments
No comments
Sign up for free and join this conversation.
If you already have a Qiita account
Why do not you register as a user and use Qiita more conveniently?
You need to log in to use this function. Qiita can be used more conveniently after logging in.
You seem to be reading articles frequently this month. Qiita can be used more conveniently after logging in.
  1. We will deliver articles that match you
    By following users and tags, you can catch up information on technical fields that you are interested in as a whole
  2. you can read useful information later efficiently
    By "stocking" the articles you like, you can search right away
ユーザーは見つかりませんでした