VBAHaskell ユーティリティを使ってみる

先日の記事VBAHaskell ユーティリティ群で書いたユーリティティモジュールは、その後いろいろ関数を追加し、記事にも反映した。
サンプルとして、このところ動きの激しい株価データにごく簡単なデータ処理をしてみよう。
まずYahoo!ファイナンスのページから、日経平均株価の時系列データをExcelシートにコピペする。

f:id:mmYYmmdd:20160130132429p:plain

このデータをもとに、前日比で大きく動いたレコードを抽出し、前日比の列を追加して表示する。
まずは配列変数へデータを取り込む。セル範囲を指定してValueプロパティを取得すれば勝手に配列になってくれる。

' データ取り込み
m = [O2:S41].Value

次に前日比を計算する。対象は「終値」の列だ。直近の日付が上になっているので、単に引き算すれば(当日値 - 前日値)となる。
行番号や列番号を指定するとき、配列のLBoundのことを気にするのは煩わしい。selectCol_b は、実際の LBound に関わりなくゼロ始まりと解釈して列を抽出する関数であり、これを使って「終値」の列番号を「4」と指定する。
adjacent_op はユーティリティ群に追加した関数で、隣り合う2項に対する演算を行うものだ。これにマイナス関数「p_minus」を渡して差分が計算できる。(当然要素数はひとつ少なくなる。)

' 前日比の列
dif = adjacent_op(p_minus, selectCol_b(m,4))

前日比の絶対値が400円を超えるものにフラグを立てよう。今作った dif 配列を「絶対値が400超」という関数でマップしてやればいい。
下の p_greater(p_abs, 400) という関数オブジェクトは見てすぐ理解できるだろう。(dumpFun 関数を使えばその構造を表示できる。)
mapF に これと dif を渡してやればフラグが作れる。

' フラグの列 : Abs(***) > 400 であれば 1、そうでなければ 0
a = mapF(p_greater(p_abs, 400), dif)

最後にやるのは、もとのデータの横に前日比の列を追加してフラグの立った行だけを抽出することだ。
前日比の列は元の行数より要素が1つ少ないので、ダミーとして"NA"という文字列を最後にくっつける。これには1次元配列(またはスカラー)を結合する関数 catV を使う。
2次元配列の列方向への結合関数 catC(CはcolumnのC)と行方向のフィルタリング関数 filterR 関数を組み合わせれば目的のものができる。

printM filterR(catC(m, catV(dif, "NA")), a)
  2016/01/29  17155.06  17638.93  16767.09   17518.3   476.849999999999
  2016/01/27  16949.19  17242.27  16947.95  17163.92   455.019999999997
  2016/01/26  16833.13  16839.52  16652.26   16708.9  -402.009999999998
  2016/01/22  16336.72  16993.96  16332.45  16958.53   941.269999999999
  2016/01/20  17030.28  17031.32  16387.61  16416.19            -632.18
  2016/01/14  17384.93  17393.83  16944.41  17240.95            -474.68
  2016/01/13  17449.12  17717.75  17414.55  17715.63   496.670000000002
  2016/01/12  17470.93  17546.57  17184.78  17218.96               -479
  2016/01/07  18139.77  18172.04  17767.34  17767.34            -423.98
  2016/01/04  18818.58  18951.12  18394.43  18450.98            -582.73
  2015/12/16   18868.2  19054.89  18859.11  19049.91   484.009999999998
  2015/12/04  19616.52   19660.9  19444.54  19504.48  -435.420000000002

浮動小数の誤差が出てカッコ悪いので、フォーマットを整えよう。

' 配列 dif を Format(, "0,00") でマップしておく
printM filterR(catC(m, catV(mapF(p_format(, "0.00"),dif), "NA")), a)
  2016/01/29  17155.06  17638.93  16767.09   17518.3   476.85
  2016/01/27  16949.19  17242.27  16947.95  17163.92   455.02
  2016/01/26  16833.13  16839.52  16652.26   16708.9  -402.01
  2016/01/22  16336.72  16993.96  16332.45  16958.53   941.27
  2016/01/20  17030.28  17031.32  16387.61  16416.19  -632.18
  2016/01/14  17384.93  17393.83  16944.41  17240.95  -474.68
  2016/01/13  17449.12  17717.75  17414.55  17715.63   496.67
  2016/01/12  17470.93  17546.57  17184.78  17218.96  -479.00
  2016/01/07  18139.77  18172.04  17767.34  17767.34  -423.98
  2016/01/04  18818.58  18951.12  18394.43  18450.98  -582.73
  2015/12/16   18868.2  19054.89  18859.11  19049.91   484.01
  2015/12/04  19616.52   19660.9  19444.54  19504.48  -435.42

まとめると、4行でできた。

m = [O2:S41].Value
dif = adjacent_op(p_minus, selectCol_b(m,4))
a = mapF(p_greater(p_abs, 400), dif)
printM filterR(catC(m, catV(mapF(p_format(, "0.00"),dif), "NA")), a)

(1/31 追記)
フォーマットを整えるためのもう少し汎用的な方法を追加した。
上の出力は小数点以下1桁の数字と2桁の数字が混在している。2桁に揃えるために、関数を列毎に適用させる columnWise_change とフォーマット指定を文字列で与える splitStr2Funs を組み合わせて以下のように書ける。(最新版の misc_utility が必要)

m = [O2:S41].Value
dif = adjacent_op(p_minus, selectCol_b(m,-1))
a = mapF(p_greater(p_abs, 400), dif)
m2 = filterR(catC(m, catV(dif, "NA")), a)
columnWise_change m2, splitStr2Funs("\\s0.00\s0.00\s0.00\s0.00\s0.00", p_str2ConvertFun, "\")
printM m2
  2016/01/29  17155.06  17638.93  16767.09  17518.30   476.85
  2016/01/27  16949.19  17242.27  16947.95  17163.92   455.02
  2016/01/26  16833.13  16839.52  16652.26  16708.90  -402.01
  2016/01/22  16336.72  16993.96  16332.45  16958.53   941.27
  2016/01/20  17030.28  17031.32  16387.61  16416.19  -632.18
  2016/01/14  17384.93  17393.83  16944.41  17240.95  -474.68
  2016/01/13  17449.12  17717.75  17414.55  17715.63   496.67
  2016/01/12  17470.93  17546.57  17184.78  17218.96  -479.00
  2016/01/07  18139.77  18172.04  17767.34  17767.34  -423.98
  2016/01/04  18818.58  18951.12  18394.43  18450.98  -582.73
  2015/12/16  18868.20  19054.89  18859.11  19049.91   484.01
  2015/12/04  19616.52  19660.90  19444.54  19504.48  -435.42

columnWise_change の第1引数は対象配列、第2引数は列ごとに適用する変換関数の配列だ。
splitStr2Funs は "\s0.00\s0.00\s0.00\s0.00\s0.00" のような文字列を区切り文字(ここでは第3引数の'\')で分けて第2引数に渡す。
第2引数の str2ConvertFun は個々の文字列("s0.00" など)を変換関数に変換する。例えば"s0.0000" だったら p_format(, "0.0000") になる。
長いので実際に使うときはまとめて1関数にしてしまえばいい。

VBAHaskell ユーティリティ群

先日の記事でVBAHaskellに数学関数を入れ始めたが、案の定その後が続かず放置状態になっている。とりあえず気にせずにユーティリティ関数群を追加する事にした。既存の関数のシンタックス・シュガーや特定の集計作業を便利にするための関数の寄せ集めだ。(適宜機能を追加していく)

シンタックス・シュガーとはいえ、3つ以上の引数を持つものを2変数にして、
mapF や zipWith、fold/scan 系の関数で使えるようになった ものもある。
モジュール名は「misc_utility」で、GithubおよびHPにある。

関数/Sub 機能 備考
p__n getNth_b(, n)の構文糖 関数オブジェクトのみ
p_typename データ型名 関数オブジェクトのみ
p_p_format Format関数 関数オブジェクトのみ
p_InStr InStr関数 関数オブジェクトのみ
p_InStrRev InStrRev関数 関数オブジェクトのみ
p_foldl1 1次元配列のfoldl1 関数オブジェクトのみ
p_foldr1 1次元配列のfoldr1 関数オブジェクトのみ
p_scanl1 1次元配列のscanl1 関数オブジェクトのみ
p_scanr1 1次元配列のscanr1 関数オブジェクトのみ
subM_R subM(m, 行範囲) の構文糖 p_subM_R を提供
subM_R_b LBound基準のsubM_R p_subM_R_b を提供
subM_C subM(m, , 列範囲) の構文糖 p_subM_C を提供
subM_C_b LBound基準のsubM_C p_subM_C_b を提供
selectRow_b LBound基準のselectRow p_selectRow_b を提供
selectCol_b LBound基準のselectCol p_selectCol_b を提供
(Sub) fillRow_b LBound基準のfillRow
fillRow_b_move LBound基準のfillRow_move
(Sub) fillCol_b LBound基準のfillCol
fillCol_b_move LBound基準のfillCol_move
- - -
adjacent_op 隣接する要素間で2項操作 p_adjacent_op を提供
(Sub) rowWise_change 2次元配列の行ごとに関数適用
rowWise_change_move 2次元配列の行ごとに関数適用 rowWise_change 後に move して返す
(Sub) columnWise_change 2次元配列の列ごとに関数適用
colomnWise_change_move 2次元配列の列ごとに関数適用 columnWise_change 後に move して返す
equal_all 全要素の等値比較 p_equal_all を提供
equal_all_pred 全要素の等値比較(述語版)
- - -
splitStr2Funs 文字列を区切って関数列へマッピング
str2SummaryFun 文字列から集計関数へ変換 group_by_partition_pointsで使用
str2ConvertFun 文字列から型変換関数へ変換
group_by_partition_points 配列のGROUP BY集計 partition_point による

p_foldl1 は、本来の foldl1 に対する関数オブジェクトではない。
foldl1は引数を3つ持つので VBAHaskell の関数オブジェクトにすることはできない。3つ目の引数はどの次元を畳み込みの軸にするかを指定している。ここでは2引数の foldl1_v という関数を実装し、それに対する関数オブジェクトとして定義している。畳み込みの軸は指定できないので1次元配列に対象が限定されるが、通常はこれで十分だと思う。
subM_RsubM_C も同様に用途を行・列に限定することで関数オブジェクトにできるようにした(p_subM_R、p_subM_C)。

adjacent_op は今までなかったタイプの関数で、隣り合うペアを関数の引数にした結果を並べる。n個の要素を持つ配列であればn-1個の要素を持つ配列が出力される。

colomnWise_change_move は2次元配列の列ごとに関数を適用して変更するもので、たとえば特定の列をString型から整数型に変換するなどに使える。
その際、適用関数の列を作らなければならないが、それをサポートするのが文字列を区切って関数列へマッピングする splitStr2Funs だ。これを使ってC言語のprintfのように"%d%d%s"といった指定が可能になる。

group_by_partition_points は GROUP BY 的な集計作業をするものだ。SQLの GROUP BY はグループを分けるカラムを決め、それ以外の集計対象カラムに対する集計方法(SUMとか)を指定する機能を持っていると思う。この関数はグループ分けをする機能は持っておらず、グループの切れ目の位置は引数で与える必要がある代わりに、全部の列が集計対象カラムになる。切れ目の位置は、Haskell_5_sortモジュールにある partition_points(またはpartition_points_pred)関数などで作ることを想定している。

' 30行6列の配列を作る(値はランダム整数)
m = makeM(30, 6, uniform_int_dist(30*6,1,8))
' 昇順ソートしておく
m = subM(m, sortIndex(m))
' 最初の列でグループ分けしたときの切れ目の位置をを求める
pp = partition_points(selectCol(m, 0))
' その切れ目でグループ分けして集計(先頭、サンプル数、SUM, Average, MIN, MAX)
g = group_by_partition_points(m, pp, "%t%c%s%a%min%max")

どういう集計をするかは3番目の変数に文字列で指定するようにした。
「% + 特定文字」が集計方法を表し、たとえば %s とか %sum が合計値といった感じだ。

集計方法 文字列
先頭の値 %t %tp %top
末尾の値 %b %btm %bottom
個数 %c %cnt %count
合計 %s %sum %計
平均 %a %avg %average
最大 %max
最少 %min

この指定には上記 splitStr2Funs を使っているが、自分で変換関数を定義することによってもカスタマイズできる。

GithubのUtilityモジュール
HPのUtilityモジュール

数学関数を入れ始める

VBAHaskellには数学関数がまだない。四則演算、指数関数、対数関数、多項式、min、max 程度が Haskell_2_stdFun にあるのと、乱数生成機が misc_random にあるだけだ。 とりあえず数学関数を少しずつ書いていくためのbasファイルを追加した。

VBA/misc_math.bas at master · mYmd/VBA · GitHub

中身はまだほとんどない。仕事で使うわけでもないからなかなか進捗しないだろう。
とりあえず Simpson法 による簡単な数値積分の関数を入れたので、そのことを書く。もちろん表面的にはループを使いたくない。

' シンプソン法による数値積分
Function integral_simpson(ByRef fun As Variant, _
                            ByVal begin_ As Double, _
                            ByVal end_ As Double, _
                            ByVal n As Long) As Double
    Dim xs As Variant, ys As Variant
    ' 独立変数の列
    xs = mapF(p_poly(, Array((end_ - begin_) / 2 / n, begin_)), iota(0, 2 * n))
    ' 関数の値を計算
    ys = mapF(fun, xs)
    ' 定数列 (1, 4, 2, 4, 2, 4, ... , 2, 4, 1) を作る
    Dim constants As Variant
    ReDim constants(0 To 2 * n)
    Call fillPattern(constants, Array(2, 4))
    constants(0) = 1
    constants(2 * n) = 1
    ' 積分値の計算
    integral_simpson = foldl1(p_plus, zipWith(p_mult, constants, ys)) * (end_ - begin_) / n / 6
End Function

4つの引数の意味は以下の通りで、返り値の型は Double だ。

引数 意味 備考
fun ByRef Variant 対象関数 VBAHaskellで定義されている関数
begin_ ByVal Douuble 積分区間の下限
end_ ByVal Douuble 積分区間の上限
n ByVal Long 分割数 それぞれの中点でも分割される

最後の備考にあるように、分割数を 100 とすると内部的には 200 分割されて、関数値は 201 個計算されることになる。Simpson法の仕様により、こうするのが自然に書ける。
念のため書いておくと、'constants' という配列は下図の「重み」に相当する部分だ。

f:id:mmYYmmdd:20151214221432p:plain

円の面積を求めてチェックする。

 ' 扇型の面積を4倍する ( 区間 [0, 1] を1,000,000 に分割して求める )
 Debug.Print 4 * integral_simpson(p_pow(p_minus(1, p_pow(,2)),0.5), 0,1,1000000)
 3.1415926534273 

 ' 真の値に近いもの
 Debug.Print 4*atn(1)
 3.14159265358979 

p_pow はべき乗を求める関数 pow_fun のファンクタで、p_pow(p_minus(1, p_pow(,2)),0.5) は √(1 - x ^ 2) を表すことになる。
100万分割で真の値と少数以下9桁まで一致した。1000分割でも4桁まで一致する。
表面的にはループはないもののスピードは遅い。100万分割だと3秒ほどかかってしまう。

中で1次元配列を他の1次元配列の繰り返しで埋める fillPattern という Sub を使っており、今回 Haskell_4_vector.bas に追加した。
ついでに多項式関数 poly もかなり遅かったので改善した。これは Haskell_2_stdFun にある。両方とも素直にループを使って実装している。このような汎用部品はパフォーマンスを優先するためしかたがない。

mbind(std::bindのメタ関数版)

std::bind は<functional>ヘッダにあるC++11の標準ライブラリで、関数やラムダ式など operator() で呼び出すことが可能なオブジェクト(Callable)に対し、引数を部分的に束縛(bind)するものだ。「呼び出し」とは実行時の関数呼び出しのことであり、引数は値である。
ここで説明する mymd::mbind はそれと同じようなことをメタ関数などのテンプレートテンプレート名に対して行うもので、型を引数に取りコンパイル時に評価される式である。

class template: mymd::mbind

namespace mymd {
    // mbindの主要インタフェース(他のメンバーは省略)
    template <template <typename...> class M, typename... binder>
    class mbind {
    public:
        template <typename... V>
            using apply = プレースホルダに型を代入してできる型;
    };
}

概要

クラステンプレート M に対し、引数を部分的に束縛(bind)する。
apply<型...> で残りの引数型を指定する。

テンプレート引数

  • M -- 複数の型を引数に取るクラステンプレート
  • binder... -- 束縛対象の型もしくはプレースホルダ(_x_, _xrr_, _xrcv_ 等)

メンバーテンプレート

  • apply -- apply<型...> としたとき型を表す式となる。ここに指定できる型の数は、mbind のテンプレート引数に指定したプレースホルダの数に等しい。すべて通常の型であった場合は M<...> の型。ひとつ以上のプレースホルダを再度指定した場合は mbind 型である。

メタ関数の例として<type_traits> の std::is_convertible をとる。

#include <type_traits>
using mymd::_x_;
using b = mymd::mbind<std::is_convertible, int, _x_>;           // 第1引数をintに束縛する
static_assert(b::apply<double>::value, "not convertible!!");    // intからdoubleへの変換は可能
static_assert(!b::apply<int*>::value, "not convertible!!");     // intからint*への変換は不可能

次はtuple の例。

#include <tuple>
// 完成型は5要素のtuple
using tuple5 = std::tuple<char, std::string, int, char, int>;
using mymd::_x_;
// 3つの型を束縛(2つはプレースホルダ)
using bind_3_5 = mymd::mbind<std::tuple, char, _x_, int, _x_, int>;
// 残り2つの型を指定する
using bind_5_5 = bind_3_5::apply<std::string, char>;
static_assert(std::is_same<bind_5_5, tuple5>::value, "!=");

定義済みプレースホルダ

applyプレースホルダに対して任意の型を代入したとき、プレースホルダの種類によって以下のように変換される。
mymd::_x_ -- 基本のプレースホルダ。型変換せずそのまま代入される。
mymd::_xrr_ -- std::remove_reference<型>::type が代入される。
mymd::_xrcv_ -- std::remove_cv<型>::type が代入される。
mymd::_xrcvr_ -- std::remove_cv<std::remove_reference<型>::type>::type が代入される。
mymd::_xdecay_ -- std::decay<型>::type が代入される。

これらはプレースホルダ生成型である mymd::_pX_ から定義されたものである。_pX_ にメタ関数の列convert...を与えて _pX_<convert...> を作れば、型 T に対してconvert0<convert1<convert2<... <T>>...> が評価されるようになっている。たとえば _xrcvr_ は _pX_<std::remove_cv, std::remove_reference>;と定義されている。

論理 or、論理 and、論理 not

mbind そのものはインスタンスを作る意味はあまりないが、インスタンスに対する operator ||, operator &&, operator ! が定義されているので、論理 or、論理 and、論理 notを作るときはインスタンス化するとやりやすい。これを型に対する式で表そうとすると相当ややこしくなってしまう。(サンプル参照)

モジュール

hppファイル: miscellaneous/mbind.hpp at master · mYmd/miscellaneous · GitHub
サンプル: miscellaneous/test_mbind.cpp at master · mYmd/miscellaneous · GitHub
[Wandbox]三へ( へ՞ਊ ՞)へ ハッハッ http://melpon.org/wandbox/permlink/6GFinX4byC5j6bCZ

P0051: C++ generic overload function

江添氏のブログに「C++標準化委員会の文書: P0050R0-P0059R0」という記事が掲載されていて、その中にオーバーロード解決のためのoverload関数テンプレートの提案が含まれている。

本の虫: C++標準化委員会の文書: P0050R0-P0059R0

オーバーロード関数を生成するための関数は自分でも作っていたが、実装をすっかり忘れているうえ名前も gen という変な関数名だったので、これを機に少し見直してみた。

(ヘッダファイル) https://github.com/mYmd/miscellaneous/blob/master/functor_overload.hpp
(サンプルcpp) https://github.com/mYmd/miscellaneous/blob/master/test_functor_overload.cpp


当初からgccでは「オーバーロードが曖昧だ」というエラーが出てどうしてもコンパイルできないという問題があり、その解消に再挑戦したのだが・・・

・・・結果、どうしてもワークアラウンドを書くことができず、再度のギブアップとなった。・・・

VBAHaskellでのコラッツ数列(その2)

以前この記事で コラッツの問題 - Wikipedia にある数列をVBAHaskellで生成するというのをやってみた。

mmyymmdd.hatenablog.com

今回は、多倍長整数を使ってもっと多くの初期値で試してみた。
このついでに以下の3つのモジュールに変更をしたので、その気があればダウンロードしてお試し願います。
VBA/misc_ratio.bas at master · mYmd/VBA · GitHub(bigIntの実装がある)
VBA/Haskell_2_stdFun.bas at master · mYmd/VBA · GitHub
VBA/Haskell_4_vector.bas at master · mYmd/VBA · GitHub

まずはコラッツ列の漸化式を定義する。
とても泥臭~いコードだが、数の大きさによって Long と 多倍長整数 bigInt をそのつど使い分けるようにしている。bigInt の計算はとても遅いからだ。

Function collatz_u(ByRef i As Variant, ByRef dummy As Variant) As Variant
    If IsNumeric(i) Then
        If 0 = i Mod 2 Then
            collatz_u = i / 2
        Else
            If i <= 715827882 Then    '715827882 = (2 ^ 31 - 2) / 3
                collatz_u = i * 3 + 1
            Else
                collatz_u = bigInt_plus(bigInt_mult(i, 3), 1)
            End If
        End If
    Else
        If bigInt_equal(bigInt_mod(i, 2), 0) Then
            collatz_u = bigInt_divide(i, 2)
            If bigInt_less(collatz_u, 2 ^ 31 - 1) Then
                collatz_u = CLng(bigInt2double(collatz_u) + 0.1)
            End If
        Else
            collatz_u = bigInt_plus(bigInt_mult(i, 3), 1)
        End If
    End If
End Function

    Function p_collatz_u(Optional ByRef firstParam As Variant, Optional ByRef secondParam As Variant) As Variant
        p_collatz_u = make_funPointer(AddressOf collatz_u, firstParam, secondParam)
    End Function

そして次に、初期値を与えてコラッツ数列を生成する関数を書く。 generate_while_notVBA/Haskell_1_Core.bas at master · mYmd/VBA · GitHub にある関数で、述語による条件が「満たされない間」繰り返し関数適用の履歴を生成するものだ。コラッツ数列の終了条件は値が 1 の項が発生することなので、述語 p_bigInt_equal(1) で表せる。

' // limit < 0 なら回数上限制約なし
Function collatz_seq(ByRef i As Variant, ByRef limit As Variant) As Variant
    collatz_seq = generate_while_not(i, p_bigInt_equal(1), p_collatz_u, limit)
End Function

    Function p_collatz_seq(Optional ByRef firstParam As Variant, Optional ByRef secondParam As Variant) As Variant
        p_collatz_seq = make_funPointer(AddressOf collatz_seq, firstParam, secondParam)
    End Function

これで出来上がりなので、とりあえず100万までやってみる。

' // 1から1000000までの整数を初期値としたコラッツ配列を生成(長時間かかる)
m = mapF(p_collatz_seq(, -1), iota(1, 1000000))
' // 列の最大の長さを求める
?foldl1(p_max, mapF(p_sizeof, m))
 525
' // 長さ525の配列を探す
?find_pred(p_equal(525), mapF(p_sizeof, m))
 837798
' // 見てみる(この記事では途中を省略)
printM m(837798)
  837799  2513398  1256699  3770098  1885049  5655148  2827574   ・・・
  ・・・ 53  160  80  40  20  10  5  16  8  4  2  1

'// m(837798)の中の最大数は何か?
?bigInt2Str(foldl1(p_bigInt_max, m(837798)))
2974984576

これは確かに VBA の Long の最大数を超えている。

qiita.com