2014年11月17日月曜日

Erlang の遅延評価で素数の無限リスト

昨日のエントリで、素数の無限リストを作る方法の基礎がわかったので、今日は Erlang と PriorityQueue を使っての試行。

Erlang の無限リストを書く場合、デフォで遅延評価ってわけじゃないので、[先頭要素|fun()-> 再帰呼出し end]として記述した上で、例えば[H|F] なら F()と明示的に呼び出して次の[H'|F']を得るという作法になるが、実装できないことはない。

これを活用して、Haskell の Data.Numbers.Primes をパクってみた。

元の Haskell コードを思い切って一言で言ってしまうと、「ホイールで生成した素数候補の無限リストを、PriorityQueue に蓄積された倍数を用いてふるいにかけつつ、同時にPriorityQueue も随時更新していくコード」といったもの。

なのでポイントとしては、ホイールと、無限リストと、PriorityQueue と、「ふるい」そのものって感じだけど、今回はホイール自体の生成は割愛し（簡単だし）ベタ書きのハードコードとした。

こんなんなった。 
-module(primes).
-compile(export_all).

enqueue(Ns, Queue) -> merge({ Ns, [] }, Queue).

merge({}, Y)  -> Y;
merge(X,  {}) -> X;
merge(X,  Y)  -> 
  case prio(X) =< prio(Y) of
    true -> join(X, Y);
    _    -> join(Y, X)
  end.

prio({[X|_], _}) -> X.
join({Ns, Qs}, Q) -> {Ns, [Q|Qs]}.

mergeAll([])  -> {};
mergeAll([Q]) -> Q;
mergeAll([Q1|[Q2|Qs]]) -> merge(merge(Q1, Q2), mergeAll(Qs)).

spin(Start, Wheel) ->
  [Start| fun() -> spin(Start, Wheel, Wheel) end].
spin(Prev, [],     Wheel) -> spin(Prev, Wheel, Wheel);
spin(Prev, [W|Ws], Wheel) ->
  Prime = Prev + W,
  [Prime | fun() -> spin(Prime, Ws, Wheel) end].

comps([N|Nl])    -> [N*N| fun()-> comps(N, Nl()) end ].
comps(P, [N|Nl]) -> [P*N| fun()-> comps(P, Nl()) end ].

sieve([N|Ns], {})    ->
  [N| fun()-> sieve(Ns(), enqueue(comps([N|Ns]), {})) end];
sieve([N|Ns], Queue) ->
  {[M|Ms], Qs} = Queue,
  if N < M -> [N| fun()-> sieve(Ns(), enqueue(comps([N|Ns]), Queue)) end];
     true  -> Q = mergeAll(Qs),
              if M == N -> sieve(Ns(),   enqueue(Ms(), Q));
                 M <  N -> sieve([N|Ns], enqueue(Ms(), Q))
              end
  end.

to_lazy([], Another) -> Another;
to_lazy([P|Ps], Another) -> [P | fun()->to_lazy(Ps, Another) end].

primes([P|Ps], Wheel) ->
  sieve(to_lazy(lists:reverse(Ps), spin(P, Wheel)), {}).
  
nth(1, [P|_]) -> P;
nth(N, [_|Pl]) -> nth(N-1, Pl()).

main(N) -> nth(N, primes([7,5,3,2], [4,2,4,2,4,6,2,6])).
実行すると以下のようになった。 
2> timer:tc(primes, main, [1000000]).
{67717516,15485863}
100万番目の素数 15485863 を得るのに、1分ちょいかかってしまっているが、まあこんなもんだ。（以前、PriorityQueue を知らずに、同じように Erlang で無限素数リストを作ってみた事があったが、それにくらべると遥かに速い。）

■ コードをざっと解説すると、まず無限リストだけど、ここでは次のものが該当する。
  • 素数候補の無限リスト: 例えば、ホイールが[2,4]だと、5から始めて[2,4,2,4,..]と足していって、[5,7,11,13,17....]を生成していく。下のコードでは、関数 spin/2, spin/3 がこれをやっている。
  • 合成数の無限リスト: 素数候補リスト中のある数と、その数以降の部分リストを用いて、合成数の無限リストを作る。例えば 7ならば、[7*7,7*11,7*13,7*17....] となる。判別された素数ごとに、この無限リストが生成される。comps/1, comps/2 がそれ。
■ 次に PriorityQueue。
Haskell の コードだけ眺めていたときは、さっぱりが訳が分からなかったけど、
sieve [3,5..] Empty から開始して 12,13個の素数を得る過程を、紙とシャーペンの筆算で評価したり図に書いたりしていると、ほぼ体得できた。

元のコードと同様、Pairing Heap 様に構成し、合成数の無限リストのそれぞれの先頭要素を比較して優先度を判定している。

型としては、元の Haskell コードでは、 
data Queue = Empty | Fork [Integer] [Queue]
といった形で定義されているけど、ここでは単純に Empty を空タプル{}、Fork を{ 倍数リスト、[倍数リスト] }として表現することにした。
enqueue/2, merge/2, prio/2, join/2, mergeAll/1 が該当する。

■ 3. そして「ふるい」。
与えられた素数候補のリストの先頭要素を、PriorityQueue の頂上にある合成数と比較して、同じなら合成数としてキャンセル、小さければ素数と判定するというコード。素数候補リスト先頭要素が、PriorityQueue 頂上の合成数より大きくなってしまう場合もあるが、その際には PriorityQueue からその合成数を除外して再トライする。

まあ、ざっくり言うとそんな感じだけど、この辺りは自然言語でコード以上に簡潔に説明するのは難しい。
投稿者 時刻: 0 件のコメント:
ラベル: ,

2014年11月14日金曜日

デンドログラムを jQuery + Underscore.js + Canvas で描いてみる

『集合知プログラミング』に載ってる、デンドログラムを書いてみる。 この ブログエントリに display:none で見えなくした div が含まれていて、『集合知プログラミング』のサンプルデータを納めてある。

データの形式は、行に「ブログ」、列に「単語」をとった表形式で、あるブログに含まれる単語の数が集計してあるというもの。

このデータをブラウザ上の Javascript で分析し、ブログ同士でピアソン相関を算出して近いものからペアにしてまとめあげていくと、クラスタができあがる。

それをツリー状に図示すると以下ようなデンドログラムになる。

ピアソン相関が 1に近づくほど（距離が小さいほど）赤みを増すように、canvas の strokeStyle を調整している。
      (function(ch03Plot01, $, undefined) {
        var $this = ch03Plot01;
        function peason(v1, v2) {
          var sums = _.chain(_.zip(v1, v2))
            .reduce(function(acc, vs){
              var a = vs[0];
              var b = vs[1];
              return [acc[0]+a, acc[1]+b, acc[2]+a*a, acc[3]+b*b, acc[4]+a*b];
            }, [0.0, 0.0, 0.0, 0.0, 0.0]).value();
          var n = v1.length;
          var num = sums[4] - sums[0]*sums[1]/n;
          var den = Math.pow((sums[2]-sums[0]*sums[0]/n)*(sums[3]-sums[1]*sums[1]/n), 0.5);

          return den==0 ? 0 : 1.0-num/den;
        }
        function Bicluster(args) {
          this.id = args.id;
          this.vec = args.vec;
          this.distance = args.distance ? args.distance : 0.0;
          this.left = args.left;
          this.right = args.right;
        }
        $this.hcluster = function (rows) {
          var distances = {};
          var currentclustid = -1;
          var clust = _.chain(_.zip(_.range(rows.length), rows))
            .map(function(idRow) {
              return new Bicluster({vec:idRow[1], id:idRow[0]});
            }).value();
          while (clust.length > 1) {
            var lowestpair = [0, 1];
            var closest = peason(clust[0].vec, clust[1].vec);
            for (var i =0; i<clust.length; ++i) {
              var x = _.range(i+10, clust.length);
              var idi = clust[i].id;
              for (var j=i+1; j<clust.length; ++j) {
                var idj = clust[j].id;
                var key = [idi, idj];
                if (!distances[key]) {
                  distances[key] = peason(clust[i].vec, clust[j].vec);
                }
                var d = distances[key];
                if (d < closest) {
                  closest=d;
                  lowestpair=[i, j];
                }
              }
            }
            var mergevec = _.chain(_.range(clust[0].vec.length))
                .map(function(i){
                  return (clust[lowestpair[0]].vec[i]+clust[lowestpair[1]].vec[i])/2.0;
                })
                .value();
            var newcluster=new Bicluster({
              vec:mergevec,
              left:clust[lowestpair[0]],
              right:clust[lowestpair[1]],
              distance:closest,
              id:currentclustid
            });
            currentclustid = currentclustid-1;
            clust.splice([lowestpair[1]], 1);
            clust.splice([lowestpair[0]], 1);
            clust.push(newcluster);
          }
          return clust[0];
        };
      } (window.ch03Plot01 = window.ch03Plot01 || {}, jQuery));
      function remove(arr, index) {
        arr.aplice(arr, index)
      }
      function getHeight(clust) {
        if (!clust.left && !clust.right) return 1;
        else return getHeight(clust.left) + getHeight(clust.right);
      }
      function getDepth(clust) {
        if (!clust.left && !clust.right) return 0;
        else Math.max(getDepth(clust.left), getDepth(clust.right))+clust.distance;
      }
      function drawDendrogram(clust, labels) {
        var h = getHeight(clust)*15;
        var w = 1200;
        var depth = getDepth(clust);
        var scalling = (w-150)/depth;

        var ctx = $('#pci-ch03-dendrogram').get(0).getContext('2d');

        $('#pci-ch03-dendrogram').attr('height', h);
        drawNode(ctx, clust, 10, h/2, scalling, labels);
        ctx.stroke();
      }

      function drawLine(ctx, x1, y1, x2, y2) {
        ctx.beginPath();
        ctx.moveTo(x1, y1);
        ctx.lineTo(x2, y2);
        ctx.stroke();
      };
      function rgbToHex(R,G,B) {return toHex(R)+toHex(G)+toHex(B)}
      function toHex(n) {
       n = parseInt(n,10);
       if (isNaN(n)) return "00";
       n = Math.max(0,Math.min(n,255));
       return "0123456789ABCDEF".charAt((n-n%16)/16)
            + "0123456789ABCDEF".charAt(n%16);
      }
      function drawNode(ctx, clust, x, y, scaling, labels) {
        if (clust.id < 0) {
          var h1 = getHeight(clust.left)*15;
          var h2 = getHeight(clust.right)*15;
          var top = y-(h1+h2)/2;
          var bottom = y+(h1+h2)/2;
          var ll = clust.distance*scaling;

          ctx.lineWidth=0.5;
          var lineColor ="#"+rgbToHex(
            Math.floor(127*(2-clust.distance*1.9)),
            Math.floor(64*(clust.distance*1.9)),
            Math.floor(127*(clust.distance*1.9)));

          ctx.strokeStyle=lineColor;

          drawLine(ctx, x, top+h1/2, x, bottom-h2/2);
          drawLine(ctx, x, top+h1/2, x+6, top+h1/2);
          drawLine(ctx, x, bottom-h2/2, x+6, bottom-h2/2);
          ctx.stroke();
          drawNode(ctx, clust.left, x+6, top+h1/2, scaling, labels);
          drawNode(ctx, clust.right, x+6, bottom-h2/2, scaling, labels);
        } else {
          ctx.fillStyle="#AAAAAA";
          ctx.fillText(labels[clust.id], (x+5), (y+2));
        }
      }
      $(function() {
        var lines = $('#blogdata').text().match(/[^\r\n]+/g);
        var colnames = _.rest(lines[0].match(/[^_]+/g));
        var $ol = $('ol');
        var size = null;
        var rowname_and_data = _.chain(_.rest(lines))
          .reduce(function(acc, line){
            var cells = line.match(/[^_]+/g);
            acc[0].push(cells[0]);
            acc[1].push(_.chain(_.rest(cells)).map(function(ch){return parseInt(ch);}).value());
            if (size == null) size = _.rest(cells).length;
            return acc;
          }, [[], []])
          .value();
        var rownames = rowname_and_data[0];
        var data=rowname_and_data[1];
        var clust = ch03Plot01.hcluster(data);
        drawDendrogram(clust, rownames);
      });

書籍では Python を使って jpg ファイルを出力していたが、ここでは計算も描画も Javascript を使ってみた。

jQuery と Underscore.js を使って、Canvas に描画している。

一昔かふた昔くらいまえは、よくこういうのに VB を使った解説とか書籍があったが、今だと Javascript + Canvas が一番お手軽かもしれない。
投稿者 時刻: 0 件のコメント:
ラベル: ,