ME差

ふとipadに入れているアプリの『血圧ノート』ってのを見たら、随分進化してた。オイラーの 自作版ソフトと同じように、平均とか偏差を出している。この結果を医者に見せてくださいって、自慢気に紹介してたぞ。

それはいいんだけど、新しい概念が発明されてた。ME差ですって。何処かの医療コンサルタントが、捻り出したんだな。薬が沢山売れて、製薬会社が儲かるようにね。

Mってのはモーニングの略、Eってのはイブニングの略。心臓がギューと縮まって血液を送り出す時、血管に圧力がかかる。これが収縮期血圧とか、最高血圧とか言われるものだ。

朝と夜でこの差が大きいと、血管プッツンして危ないですよってのが、ME差らしい。大体、この差が15から20ぐらいより大きいと、問題ですよって脅かす訳だ。幅が有るのは、まだ提唱中だから、世間一般に認知されていないからだろうね。

これと、朝晩の血圧の平均値の135を組み合わせる。135ってのは、高血圧学会なりWHOなりがお墨付きを与えた、高血圧の指標だ。

           | 平均135以下  | 平均135越え
-----------+-----------------+----------------
ME差15越え | 早朝血圧上昇型 | 早朝優位高血圧
ME差15以下 | 正常血圧     | 持続性高血圧

起きる時はゆっくりと、気温差を無くす、早朝高血圧に適した薬の服用(これを言いたいのが見え見え)。某タレントさんが罹患した、悪性高血圧は、、、と、脅かすのに余念がないな。

そういう発明に乗って、オイラーも調べてみるか。幸いたっぷりとデータが有るからね。 時系列にため込んだCSVファイルになってる。テキストなんで、awkと相性が良いでしょう。 ってんで、awkを思い出しながら、スクリプトを書いてみた。

旅行に行ってるとかして、朝晩のデータが欠落してる場合が有るので、それは除外するようにしたよ。ってか、そうしないといけない。まあ、朝晩のサンプリング。かっこよく言うと、離散的だな。

腹筋とかやると血圧が上昇する事を、体で感じるんだけど、こういうのはどうなのよ。血圧上昇するから、そんな無茶するなって事?

#! /usr/bin/awk -f
BEGIN { FS = ","; }

{
    hh = substr($1,7,2);
    ymd= substr($1,1,6);
    if ( hh < 12 ){       # Case AM
        mb = $2;
        today = ymd;
        next;
    }
    if (today == ymd){    # Case PM
        eb = $2;
        av = (mb + eb) / 2;
        dl = mb - eb;
        if (av >  135 && dl > 15)   r1 +=1;
        if (av <= 135 && dl > 15)   r2 +=1;
        if (av >  135 && dl<= 15)   r4 +=1;
        if (av <= 135 && dl<= 15)   r3 +=1;
        ttl++;
    }
}

END {
    printf "%8d%8d\n%8d%8d\n", r2,r1,r3,r4;
    printf " Total count = %d\n", ttl;
}

実行時間を測りつつ、走らせてみた。

[ham@cent tmp]$ time ./me.awk current.csv
     856       4
    1773      75
 Total count = 2708

real    0m0.012s
user    0m0.010s
sys     0m0.001s

リナ陣営のやつは便利な機能が搭載されてるんだけど、知らず知らずのうちにそれを使ってしまって動かない事があるので、shだとかawkのコード書きはBSDで行うようにしてる。

by octave

丁度、octaveにハマっているんで、同じ事をoctaveでやってみる。

% me blood
tic();
bld = csvread('current.csv');
am  = bld( rem(bld(:,1),100)<12, [1 2] ); % [ [ymdh hi lo ps] ... ]
pm  = bld( rem(bld(:,1),100)>12, [1 2] );
am  = [fix(am(:,1)/100) am(:,2)];         % [ [ymd hi] ... ]
pm  = [fix(pm(:,1)/100) pm(:,2)];
tbl = pm(:,1);                            % for lookup (must be sorted)
toc()
res = zeros(2,2);
for i=1:length(am)
   if ( x = lookup(tbl, am(i,1), 'm') )
     av = (am(i,2) + pm(x,2)) / 2.0;
     dl = am(i,2) - pm(x,2);
     if (av >  135 && dl > 15)   res(1,2) +=1; endif
     if (av <= 135 && dl > 15)   res(1,1) +=1; endif
     if (av >  135 && dl<= 15)   res(2,2) +=1; endif
     if (av <= 135 && dl<= 15)   res(2,1) +=1; endif
   endif
endfor
res
sum(sum(res))
toc()

tic,tocを忍ばせて、経過時間をざっくりとモニターするようにした。

[ham@cent tmp]$ octave me.m
warning: suggest parenthesis around assignment used as truth value near line 12, column 11 in file '/tmp/me.m'
Elapsed time is 0.00973105 seconds.
res =

    856      4
   1775     75

ans =  2710
Elapsed time is 0.127146 seconds.

結果がawkのそれと異なっているけど何故? 気になるぞ! こういう動くけど、結果が違うってバグを見つけるのは至難の技が必要そう。 別の方法を試してみるかな。

その前に、今回の作の意図をば。なるべく、らしく書くのが目標だった。 前半は、らしく書けたけど、後半はありきたりだなあ。そして、圧倒的に後半が遅い! forを使ってるから? いいや、無理してlookupなんてのを使ったからかな。

octave:2> size(am)
ans =

   2723      2

octave:3> size(pm)
ans =

   2718      2

朝と晩の測定個数が違う。両方のデータが揃っているのを対象にしたかったんで、晩の日付のテーブルを作り、その中に朝の日付が含まれているかってのを一々lookupしてる。見つかれば、そのindexが得られる。見つからないと0が返る。こういう手法で、同時に存在するデータを得てる。超無駄な事だ。

前半部分も、手軽に何個もマトリックスを作っているけど、場所の確保やら値の設定で無駄を してるな。手軽に出来るからこそ、戒めないとな。

csvread()にどれぐらいの時間がかかっているのだろう?csvread()の後ろにtoc()を置いてみた。

octave> me
Elapsed time is 0.00999403 seconds.
Elapsed time is 0.01033 seconds.

前半部分で、ほとんどの時間はcsvread()に費やされているとな。IOは高価なオペレーションなんだな。マトリックスを対象にした演算は気にする程でないとな。

まてまて、遅いのはIOがからんでいるって、そんな大雑把な事で委員会? csvreadって、読み取った数字文字列をバイナリー表現の数値に変換してからマトリックスに格納するんだよな。 何行あるか分からんファイルを相手にするんだから、一行分づつをappendするんだかLispっぽくconsするんだかは知らないけど、高負荷な事をやってるはずだな。

csvread()

何か高速にする方法が無いかと、help csvread する。dlmreadへのラッパーって事が判明。

m/io/csvread.m

function x = csvread (filename, varargin)
  x = dlmread (filename, ",", varargin{:});
endfunction

続いて、help dlmread する。指定した所だけ読み込む方法が解説されてた。いわゆる範囲を指定出来る。それもあのEXCELみたいにね。そんな事で、スピード比べ。

clear all
tic()
% bld = csvread('current.csv');
bld = dlmread('current.csv', ',' , 'A1..B999999');
toc()

A列の1行目からB列の99999行分読み取ってくれ。99999ってのは十分と言う意味ね。これで、使いもしない列はオミットされるから、速くなる期待。

octave> Elapsed time is 0.031836 seconds.   % by csvread
octave> Elapsed time is 0.032347 seconds.
octave> Elapsed time is 0.030863 seconds.
octave> Elapsed time is 0.0301418 seconds.  % by dlmread
octave> Elapsed time is 0.02987 seconds.
octave> Elapsed time is 0.0301368 seconds.

3回づつ、測定してみた。

期待したんだけど、誤差の範囲内だな。データはint型で格納されてると思ったんだけど、double型だった。デフォルトの数値型はdoubleになってる。今まで気がつかなかったぞ。

octave:40> s="19070121,123\n888,56789\n"
s = 19070121,123
888,56789

octave:41> q = textscan(s, "%d,%d");
octave:42> q
q =
{
  [1,1] =

    19070121
         888

  [1,2] =

      123
    56789

}

こういう機能が有るから使え? 大カッコで括られていて、cellクラスとか言うらしい。深入りはしないでおこう。

検証

で、一番の問題は、正しい答えが出てこない事。どんなにスピードが速くても、間違っていたんじゃ、目もあてられない。

全てを疑って、元データから検証してみる。朝晩に測定値の有る日数の確認。cutコマンドのバイト位置による抽出でymdを取り出し、uniq して、その個数を数える。

[ham@cent tmp]$ cut -b1-6  current.csv | sort |  uniq -c | grep ' 2 ' | wc
   2706    5412   40590
[ham@cent tmp]$ cut -b1-6  current.csv | sort |  uniq -c | sort
         :
      2 190630
      3 131013
      3 131021

自信のある(?)awkの実行結果とも2日程ズレていたので、もしやと思って眺めてみると、だぶりを発見。こういうのに関しては、自作のアプリでは検出出来ないからなあ。

[ham@cent tmp]$ grep 131013 current.csv
13101303,116,63,56
13101304,106,64,59
13101321,107,56,63
[ham@cent tmp]$ grep 131021 current.csv
13102103,116,72,55
13102104,116,73,58
13102121,126,70,60

入力を間違ったのかな? 修正しておこう。

この修正を施したら、awkが出す結果とoctaveの出す結果が一致した。第三者の眼って大事だなあ。大事に至らなくてよかったわい。

profile

プロファイラーが用意されてるので、そいつのお世話にもなってみる。

octave:1> profile on
octave:2> me
octave:3> profile off
octave:4> profshow
   #  Function Attr     Time (s)   Time (%)        Calls
--------------------------------------------------------
   1        me             0.112      83.99            1
   3   dlmread             0.010       7.16            1
  11    lookup             0.006       4.67         2721
   6  binary >             0.002       1.27         8125
  14 binary <=             0.001       1.09         8124
  12  binary +             0.001       0.64         2708
  13  binary -             0.001       0.49         2708
   7  binary /             0.001       0.40         2710
  16      disp             0.000       0.11            2
   4       rem             0.000       0.07            2
   2   csvread             0.000       0.04            1
  18   profile             0.000       0.03            1
  15   display             0.000       0.01            2
   8       fix             0.000       0.01            2
   5  binary <             0.000       0.01            1
  17       sum             0.000       0.01            2
   9     zeros             0.000       0.00            1
  19    nargin             0.000       0.00            1
  21     false             0.000       0.00            1
  20 binary !=             0.000       0.00            1

profile('clear')で、ログデータがクリア出来る。一定時間毎に割り込んで、何をしてたか 調べているんだろうから、数度実行してみると良いかも知れない。

octave:6> profexplore

Top
  1) me: 1 calls, 0.134 total, 0.112 self
  2) profile: 1 calls, 0.000 total, 0.000 self

profexplore> 1

Top
  me: 1 calls, 0.134 total, 0.112 self
    1) csvread: 1 calls, 0.010 total, 0.000 self
    2) lookup: 2721 calls, 0.006 total, 0.006 self
    3) binary >: 8125 calls, 0.002 total, 0.002 self
    4) binary <=: 8124 calls, 0.001 total, 0.001 self
    5) binary +: 2708 calls, 0.001 total, 0.001 self
    6) binary -: 2708 calls, 0.001 total, 0.001 self
    7) binary /: 2710 calls, 0.001 total, 0.001 self
    8) display: 2 calls, 0.000 total, 0.000 self
    9) rem: 2 calls, 0.000 total, 0.000 self
    10) fix: 2 calls, 0.000 total, 0.000 self
    11) binary <: 1 calls, 0.000 total, 0.000 self
    12) sum: 2 calls, 0.000 total, 0.000 self
    13) zeros: 1 calls, 0.000 total, 0.000 self
    14) length: 1 calls, 0.000 total, 0.000 self

こちらは、ターゲットを絞った内容確認になる。

ここには出ていないけど、forで回すのが決定的に遅い。forを何とか駆逐出来ないか?

for less

強引にforを表舞台から消してみた。 その仕掛けは、arrayfunと言う関数。アレー(複数指定可)と、各エレメントに対して実行したい関数のハンドルを指定する。

gaucheあたりのforeachだな。関数のハンドルってoctaveの方言が出てきたけど、scheme風に言うと、lambdaだ。

下記で定義した関数の引数xに対して処理を書いていく。関数なので、基本的には関数外は見えない。それじゃ困ると言う事で、関数の内外で、グローバル宣言と言う怪しげ事をする。超不格好で、見るからにこんな事をするんじゃ無いってのがプンプンする。

朝と晩で共通に表れる日を、集合演算 intersectで求めてみた。こんな事も出来るよって見本かな。

clear -a
bld = csvread('current.csv');
am  = bld( rem(bld(:,1),100)<12, [1 2] );   % [ ymdh hi lo ps; ... ]
pm  = bld( rem(bld(:,1),100)>12, [1 2] );
global am2  = [fix(am(:,1)/100) am(:,2)];   % [ ymd hi; ... ]
global pm2  = [fix(pm(:,1)/100) pm(:,2)];
global tba = am2(:,1);                      % [ ymd ... ]
global tbp = pm2(:,1);
both= intersect(tba,tbp);

global res = zeros(2,2);

function diside (x)
  global am2 pm2 tba tbp res
  i = lookup(tba, x, 'm');  mb = am2(i, 2);
  j = lookup(tbp, x, 'm');  eb = pm2(j, 2);
  av = (mb + eb) / 2;
  dl = mb - eb;
  if (av >  135 && dl > 15)   res(1,2) +=1; endif
  if (av <= 135 && dl > 15)   res(1,1) +=1; endif
  if (av >  135 && dl<= 15)   res(2,2) +=1; endif
  if (av <= 135 && dl<= 15)   res(2,1) +=1; endif
endfunction

arrayfun(@diside, both)
res
sum(sum(res))

気になる実行時間だけど、最初に書いたものより遅くなった。Elapsed time is 0.147579 seconds. ですってさ。多分、lookupが足を引っ張っているんだろうね。

octave:2> whos
Variables in the current scope:

   Attr Name        Size                     Bytes  Class
   ==== ====        ====                     =====  =====
        am       2721x2                      43536  double
     g  am2      2721x2                      43536  double
        ans         1x1                          8  double
        bld      5439x4                     174048  double
        both     2708x1                      21664  double
        pm       2718x2                      43488  double
     g  pm2      2718x2                      43488  double
     g  res         2x2                         32  double
     g  tba      2721x1                      21768  double
     g  tbp      2718x1                      21744  double

Total is 51664 elements using 413312 bytes

スクリプトの冒頭に入れている clear -a は、グローバル宣言した変数も含めて、クリアする為だ。こいつを入れておかないと、結果が累積されて、ホェーとなる。

res =

    5136      24
   10638     450

ans =  16248
global res = zeros(2,2);
global today mb

function diside (x,y)
  global today mb res
  ymd = int32(x) / 100;
  if (rem(x,100) < 12)
    today = ymd;
    mb = y;
  elseif (today == ymd)
    av =( mb + y) /2;
    dl = mb - y;
     if (av >  135 && dl > 15)   res(1,2) +=1; endif
     if (av <= 135 && dl > 15)   res(1,1) +=1; endif
     if (av >  135 && dl<= 15)   res(2,2) +=1; endif
     if (av <= 135 && dl<= 15)   res(2,1) +=1; endif
  endif
endfunction

bld = csvread('current.csv');
arrayfun(@diside, bld(:,1), bld(:,2))
res
sum(sum(res))

またawkのそれに似せて書いたスクリプトだと、 Elapsed time is 0.165908 seconds. こんな具合に更に遅くなった。隠れforは元凶だな。

forが遅いのか? forループの中でインタプリットしてプリミティブ関数を呼び出すのが遅いのか? それが問題だ。

roots

ME差を肴に、色々やったけど、ちょっと趣向を変えて、4つのエリアに分かれるその原点を満たす血圧はどんな値になるんだろう? mbを朝血圧、ebを夜血圧とすると

mb + eb = 270
mb - eb = 15

説明では、平均って事で、2で割ってたけど、それを払ってある。こうすると、連立方程式になるな。

左辺の式において、係数だけを抜き出した行列aを作る。右辺の値を列ベクトルで表す。そうしておけは、いろいろな方法で、連立方程式を解けるぞ。(一つ条件がある。det(a) != 0 の場合だけだよ)

octave:1> a  = [1 1; 1 -1]
a =

   1   1
   1  -1

octave:2> b = [270; 15]
b =

   270
    15

octave:3> inv(a) * b
ans =

   142.50
   127.50

octave:6> linsolve (a,b)
ans =

   142.50
   127.50

octave:7> a \ b
ans =

   142.50
   127.50

N次式の解法も出来る。helpよりの抜粋。

          p(x) = x^2 - 5.

          c = [1, 0, -5];
          roots (c)
          =>  2.2361
          => -2.2361

別な例

octave:10> roots([2 3 -4 5])
ans =

  -2.62482 + 0.00000i
   0.56241 + 0.79759i
   0.56241 - 0.79759i

2*x^3 + 3*x^2 - 4*x + 5 = 0 の根を求めるんだな。便利、って言っても、身の回りで こんなのの出番は有るのか? (複素数なんて普通の人は無用の長物だ)