レインフロー法 #2（デバッグ準備編）

レインフロー法 #1（準備編） - potass' blog の続きで Python に
レインフローモジュールがあったため、次回の実装編でデバッグが出来るように触ってみる。
加えて実装上のアルゴリズムがどうなっているかも取り上げる。
pypi.org

インストール

$ pip install rainflow

で無事インストール（Python 3.8.6、Windows10 64bit）。

アルゴリズム

このモジュールの使い方にもあるようにcount_cycles(series, ndigits, nbins, binsize)（各パラメタについては下記とのこと）

Parameters
    ----------
    series : iterable sequence of numbers
    ndigits : int, optional
        Round cycle magnitudes to the given number of digits before counting.
        Use a negative value to round to tens, hundreds, etc.
    nbins : int, optional
        Specifies the number of cycle-counting bins.
    binsize : int, optional
        Specifies the width of each cycle-counting bin

    Arguments ndigits, nbins and binsize are mutually exclusive.

    Returns
    -------
    A sorted list containing pairs of range and cycle count.
    The counts may not be whole numbers because the rainflow counting
    algorithm may produce half-cycles. If binning is used then ranges
    correspond to the right (high) edge of a bin.

が入出力の際にメインとなる関数であるが実態（計算）はそこで呼ばれているextract_cycles(series)が基本となる。
さて、extract_cycles(series)では

与えられたデータ列seriesをreversal(series)によって極値のみのデータ列に変換。*1 変換方法は隣接するデータ列の差の掛け算が負になる場合はその点を抽出、そうでない場合は無視するというシンプルなもの。
先の極値のみのデータ列に対し先頭のデータから順番に読み込んでいき、下図の2パターンになったときは該当するやり方で小ループを除去する。

f:id:potass:20210103070615j:plain:h300
が実行されており、前回記載した P/V 差法そのものである。

お試し

お試しコードは下記。実行自体はrainflow.extract_cycles(sgnl)のみ。import rainflowをお忘れなく。

gist.github.com

実行結果は下記。ピークは scipy で抽出したもの。これがrainflow.extract_cycles(sgnl)の出力と対応していることを確認できた。
次回は Excel VBA での実装。

f:id:potass:20210103084849j:plain:h300

[[ 0.20100503  0.92686369]
 [ 1.00502513  1.97208224]
 [ 1.76884422  2.34561486]
 [ 2.53266332  1.40625694]
 [ 3.33668342  0.76749291]
 [ 0.56281407 -1.01860648]
 [ 1.36683417 -0.0552989 ]
 [ 2.17085427 -0.44787773]
 [ 2.95477387 -1.48082349]
 [ 3.71859296 -1.45553457]]
-----------------------------------------
range, mean, count, i_start, i_end
1.326863691183818 0.26343184559190913 0.5 0 10
1.945470171916998 -0.04587139477468083 0.5 10 28
2.0273811415975134 0.9583916687347154 1.0 50 68
1.854134671396499 0.479189607932012 1.0 108 126
3.364221339924992 0.663504189229316 0.5 28 88
3.826438349059406 0.43239568466210876 0.5 88 147
2.2483163957095447 -0.35666529201282193 0.5 147 166
2.2230274803187586 -0.3440208343174289 0.5 166 185
2.101825614128834 -0.40462176741239136 0.5 185 199

*1:この変換で吐き出されるインデックスの番号が1つ分ズレてる。どうせその後使用してないから別にいいのだが。ただ、enumerate(series, start=1) としていて意図的にも感じるので何か意味があるのかもしれない…。