似たようなことがしてみたい人へ

えー、一応。似たようなことをしたい人がいる可能性を考えて、参考資料を。まずは上司/スポンサー説得資料。

基本的にこれらの本に書いてあるのとやっていることは同じである、と主張しても大きく嘘にはならないと思います。いや、本当に同じ様なことをしているんですけれどね。統計学を仕事に応用しているんだし。ただ、あまりすさまじくかっこいい具合にありえないようなものを発見できる、と言う風に誤解されるとそれはそれで後が大変でしょうから、こう…加減と言うものは各自で調整してください。

で、これで何を手に入れるかと言うと、もちろん「計算機」です。大事なのは「64bit」と「メモリがたくさん」。いまどきのx86マシンならば EM64T対応な段階で SSEとかも相応に早いです。もちろん、Core2Duo とか Core2Quad とか Corei7 とかがいいわけですが、動作クロックは最低でも構わないかと。むしろ大事なのはメモリ。

Windows はやめておいた方がよいです。メモリの利用効率が悪い、あれは。ファイルキャッシュのデザインが unix 以前(1969年以前、と言うこと)のOSそっくりだ。計算用データを保持するために1pageでも余計にメモリを使いたいときに、ガッツリ kernel 内部のファイルキャッシュがページを握って離さない、というWindows NT/Longhorn デザインは勘弁して欲しい(つーか Windows7 では是非その辺りも改善していただきたいものだ)。と言うわけで、とりあえずお勧めは Linux 。マルチコア、HT対応、Intel & AMD 御謹製のプロセッサコントロールコードが使えるのは Linux ぐらいですから。

次に。あなた自身のための入門本。

多分、この一冊と Wikipedia 辺りが入門には丁度良いのではないかと。もちろん、前出の2冊も読んでいただくのですがあれは「どういう場面で統計を使うか」なのに対して、ここからは「その場面で何を使うか」。

で、最後に「なんとなく判ったけど、コーディングが面倒だな」とか「もっと詳しく知りたいが、できれば数学的なほうじゃなくて…」という人のための本。

多分これぐらいあれば、1年はここに書いてある内容で遊べます。

えぇ、「遊んで暮らせます」というオチだったらよかったんですけどねぇ。残念ながらそこまでは…

リファクタリング-5-

さて。一番重たい部分を計算するプログラム公開の前に。どうやら問題を理解してもらえていないようなので、もう一度書いておく。

m 個のデータ列がある。各データ列は N 個のデータからなる。N個のデータはm個のデータ列の間で互いにどれとどれが対応するのかは判っているとする。

これらのデータ列同士の相関を「総当りで」求めたい。相関係数 r は r(x,y) と r(y,x) （ただし、x, y はデータ列) は同じ値になるので、m個のデータ列同士の総当たり戦は m*(m-1)/2 通りになる。これを

面倒なコーディングを可能な限り行わない

形で、実用的な速度で解きたい。

R であれ、ライブラリであれ、モジュールであれ、1組2個のデータ列同士の相関を求める関数は存在するし、それはどれを使っても相応に早い。しかし「総当たり戦」の場合、その関数を m*(m-1)/2 回呼び出さなくてはいけない。

相関を求める式には、データ列単位で計算しておけばそれを使いまわせるものが結構ある。平均値を求めたり、各データと当該平均値との差分を求めたり、それら差分の二乗値の和の平方根を取ったりした結果などがそうだ。

ある任意のデータ列は、(m-1)/2 回相関を求めるために呼び出される。と言うことは、上記の値は、最初に1回だけ計算しておけば ((m-1)/2) - 1 回分、計算をしなくてよい。

このような最適化を施したライブラリは、いまだに存在しない。というか、存在しなくは無いのだが、メモリを食いすぎるため、私のPCでは遅すぎるのだ(32bit環境だと、そもそもアドレス空間が足りなくて core 吐いて終わる。64bit 環境は、メモリが足りないため swap IO ばかりで計算が全然前に進まない)。計算環境が64bit環境だけで作れて、メモリが潤沢にあれば問題は無いのだが…

その一方で。この相関は何のためにとるのかと言うと、「Windows の perfmon の出力から、何か性能劣化の兆候が見えないか?」とか「性能が劣化しているのだが、何が原因かわからないか?」という分析をする際の、調査対象に優先度をつけるために行っている。つまり計算誤差なんぞ多少あっても構わないのだ。まじめな数値計算ではなく、誤差が多少あってもいいから高速に結果を出して欲しい。

実用になる範囲であれば、perl 辺りでライブラリを叩けば十分だったのだが、速度が遅すぎる以上 C でコーディングをするしかないわけだ。

ちなみに、mは 27000．Nが 675。 perl で組んでいたときは、1組当たり 9msec かかっていた。つまり37.9673.. 38日かかる計算になる。これが C だと 37分…まぁIOの関係で40分ぐらい。「値が変化しない」とわかりきっている部分が半分近くあり、それらについては相関は絶対 NaN になると判っているのでメモリも消費しなくなった。すると32bit環境でも一発で計算できるではないか…。NIH(Not Invented Here)症候群と言われようが、自分でコーディングする気になるというものだよ。

次は一番重たい部分を計算するCのプログラム

correlation.c

#include <stdio.h>
#include <stdlib.h>
#include <string.h>
#include <unistd.h>

#include <math.h>

#define        DEBUG   (0)



typedef struct __data {
 struct __data        *next;
 struct __data        *prev;
 char                *element_name;
 int                        correlation;
 double                sum;
 double                mean;
 double                sqsum;
 double                subpart;
 long                numdata;
 double                *datalist;
} data;


data*        read_each( char *listfilename );
data*        read_datafile( char *linebuffer );
void    add_data_array( data **datalistpp, data *new_datap );
void    calculate( data *data_list );
void        remove_trailing_newline( char *linebuffer );

void        help( char *programname );
void        *o_malloc( size_t size, const char *msg, const char *func, const long line );



int
main( int argc, char *argv[] )
{
 data        *data_list;
 if ( argc != 2 ) {
     help( argv[0] );
 }

 data_list        = read_each( argv[1] );
 if ( data_list == NULL ) {
     fprintf( stderr, "no data given.\n" );
     exit( 1 );
 }

 calculate( data_list );
}

data*
read_each( char *listfilename )
{
 FILE        *fp;
 char        *linebuffer;
 int                linebufsize;
 data        *read_datap        = NULL;


 linebufsize        = getpagesize();
 linebuffer        = o_malloc( linebufsize, "unable to prepare linebuffer\n", __func__, __LINE__ );

 fp        = fopen( listfilename, "r" );
 if ( fp == NULL ) {
     fprintf( stderr, "unable to open file %s\n", listfilename );
     exit( 1 );
 }

 while( fgets( linebuffer, linebufsize, fp ) != NULL ) {
     data        *new_datap;
     long        len;

     if ( linebuffer[0] == '#' ) {
         continue;
     }

     remove_trailing_newline( linebuffer );

     new_datap        = read_datafile( linebuffer );
     if ( new_datap->correlation != 0 ) {
         add_data_array( &read_datap, new_datap );
     } else {
         /* take it away from the memory.*/
         free( new_datap->datalist );
         free( new_datap->element_name );
         free( new_datap );
     }
 }
 fclose( fp );

 return read_datap;
}

data*
read_datafile( char *datafilename )
{
 data        *newdata;
 FILE        *fp;
 char        *linebuffer;
 size_t        linebufsize;
 char        *tagname;


 newdata        = o_malloc( sizeof( data ), "unable to allocate memory for newdata\n", __func__, __LINE__ );

 linebufsize        = getpagesize();
 linebuffer        = o_malloc( linebufsize, "unable to prepare linebuffer for read_datafile()\n", __func__, __LINE__ );

 fp        = fopen( datafilename, "r" );
 if ( fp == NULL ) {
     fprintf( stderr, "unable to open file %s as data file target\n", datafilename );
     exit( 1 );
 }


 /* read tagname */
 {
     char        *p;
     long        dataoffset;

     if( fgets( linebuffer, linebufsize, fp ) == NULL ) {
         fprintf( stderr, "unable to read tagname at file %s\n", datafilename );
         exit( 1 );
     }
     remove_trailing_newline( linebuffer );

     p        = strstr( linebuffer, "tagname\t" );
     if ( p != linebuffer ) {
         fprintf( stderr, "unable to find \"tagname\t\" indicator in %s\n", datafilename );
         exit( 1 );
     }
  
     dataoffset        = strlen( "tagname\t" );

     newdata->element_name        = strdup( linebuffer + dataoffset );
     if ( newdata->element_name == NULL ) {
         fprintf( stderr, "unable to copy tagname %s in %s\n",
                  linebuffer + dataoffset, datafilename );
         exit( 1 );
     }
#if DEBUG
     fprintf( stderr, "tagname was %s\n", newdata->element_name );
#endif/*DEBUG*/
 }

 /* read Correlation */
 {
     char        *p;
     long        dataoffset;

     if( fgets( linebuffer, linebufsize, fp ) == NULL ) {
         fprintf( stderr, "unable to read Correlation at file %s\n", datafilename );
         exit( 1 );
     }
     p        = strstr( linebuffer, "Correlation\t" );
     if ( p != linebuffer ) {
         fprintf( stderr, "unable to find \"Correlation\t\" indicator in %s\n", datafilename );
         exit( 1 );
     }
     dataoffset        = strlen( "Correlation\t" );
     if ( strncasecmp( linebuffer + dataoffset, "OK", 2 ) != 0 ) {
         newdata->correlation        = 0;
     } else {
         newdata->correlation        = 1;
     }
#if DEBUG
     fprintf( stderr, "correlation was %d\n", newdata->correlation );
#endif/*DEBUG*/
 }


 /* read sum */
 {
     char        *p;
     long        dataoffset;
     int        res;

     if( fgets( linebuffer, linebufsize, fp ) == NULL ) {
         fprintf( stderr, "unable to read sum at file %s\n", datafilename );
         exit( 1 );
     }
     p        = strstr( linebuffer, "sum\t" );
     if ( p != linebuffer ) {
         fprintf( stderr, "unable to find \"sum\t\" indicator in %s\n", datafilename );
         exit( 1 );
     }
     dataoffset        = strlen( "sum\t" );
     newdata->sum        = 0.0;
     res        = sscanf( linebuffer + dataoffset, "%lf", &( newdata->sum ));
     if ( res != 1 ) {
         perror( "sscanf failed when reading sum.\n" );
         exit( 1 );
     }
#if DEBUG
     fprintf( stderr, "sum was %lf\n", newdata->sum );
#endif/*DEBUG*/
 }


 /* read mean */
 {
     char        *p;
     long        dataoffset;
     int        res;

     if( fgets( linebuffer, linebufsize, fp ) == NULL ) {
         fprintf( stderr, "unable to read mean at file %s\n", datafilename );
         exit( 1 );
     }
     p        = strstr( linebuffer, "mean\t" );
     if ( p != linebuffer ) {
         fprintf( stderr, "unable to find \"mean\t\" indicator in %s\n", datafilename );
         exit( 1 );
     }
     dataoffset        = strlen( "mean\t" );
     newdata->mean        = 0.0;
     res        = sscanf( linebuffer + dataoffset, "%lf", &( newdata->mean ));
     if ( res != 1 ) {
         perror( "sscanf failed when reading mean.\n" );
         exit( 1 );
     }
#if DEBUG
     fprintf( stderr, "mean was %lf\n", newdata->mean );
#endif/*DEBUG*/
 }


 /* read sqsum */
 {
     char        *p;
     long        dataoffset;
     int        res;

     if( fgets( linebuffer, linebufsize, fp ) == NULL ) {
         fprintf( stderr, "unable to read sqsum at file %s\n", datafilename );
         exit( 1 );
     }
     p        = strstr( linebuffer, "sqsum\t" );
     if ( p != linebuffer ) {
         fprintf( stderr, "unable to find \"sqsum\t\" indicator in %s\n", datafilename );
         exit( 1 );
     }
     dataoffset        = strlen( "sqsum\t" );
     newdata->sqsum        = 0.0;
     res        = sscanf( linebuffer + dataoffset, "%lf", &( newdata->sqsum ));
     if ( res != 1 ) {
         perror( "sscanf failed when reading sqsum.\n" );
         exit( 1 );
     }
#if DEBUG
     fprintf( stderr, "sqsum was %lf\n", newdata->sqsum );
#endif/*DEBUG*/
 }


 /* read subpart */
 {
     char        *p;
     long        dataoffset;
     int        res;

     if( fgets( linebuffer, linebufsize, fp ) == NULL ) {
         fprintf( stderr, "unable to read subpart at file %s\n", datafilename );
         exit( 1 );
     }
     p        = strstr( linebuffer, "subpart\t" );
     if ( p != linebuffer ) {
         fprintf( stderr, "unable to find \"subpart\t\" indicator in %s\n", datafilename );
         exit( 1 );
     }
     dataoffset        = strlen( "subpart\t" );
     newdata->subpart        = 0.0;
     res        = sscanf( linebuffer + dataoffset, "%lf", &( newdata->subpart ));
     if ( res != 1 ) {
         perror( "sscanf failed when reading subpart.\n" );
         exit( 1 );
     }
#if DEBUG
     fprintf( stderr, "subpart was %lf\n", newdata->subpart );
#endif/*DEBUG*/
 }


 /* read datas */
 {
     char        *p;
     long        dataoffset;
     int        res;

     if( fgets( linebuffer, linebufsize, fp ) == NULL ) {
         fprintf( stderr, "unable to read datas at file %s\n", datafilename );
         exit( 1 );
     }
     p        = strstr( linebuffer, "datas\t" );
     if ( p != linebuffer ) {
         fprintf( stderr, "unable to find \"datas\t\" indicator in %s\n", datafilename );
         exit( 1 );
     }
     dataoffset        = strlen( "datas\t" );
     newdata->numdata        = 0;
     res        = sscanf( linebuffer + dataoffset, "%ld", &( newdata->numdata ));
     if ( res != 1 ) {
         perror( "sscanf failed when reading datas.\n" );
         exit( 1 );
     }
#if DEBUG
     fprintf( stderr, "datas was %d\n", newdata->numdata );
#endif/*DEBUG*/
 }

 /* read data list */
 {
     long i;
     int        res;

     newdata->datalist
         = malloc( sizeof( double ) * newdata->numdata );
     if ( newdata->datalist == NULL ) {
         fprintf( stderr, "not enough memory for datalist: %s\n", datafilename );
         exit( 1 );
     }

     if( fgets( linebuffer, linebufsize, fp ) == NULL ) {
         fprintf( stderr, "unable to read empty line at file %s\n", datafilename );
         exit( 1 );
     }

     for ( i = 0; i < newdata->numdata; i++ ) {
         if( fgets( linebuffer, linebufsize, fp ) == NULL ) {
             fprintf( stderr, "unable to read datas at file %s\n", datafilename );
             exit( 1 );
         }
         res        = sscanf( linebuffer, "%lf", &( newdata->datalist[i] ));
         /* if we failed in reading data, that means that element should be 0.0 */
         if ( res != 1 ) {
             newdata->datalist[i]        = 0.0;
         }

#if DEBUG
         fprintf( stderr, "datas was %lf\n", newdata->datalist[i] );
#endif/*DEBUG*/
     }
 }

 fclose( fp );
 return newdata;
}

void
add_data_array( data **datalistpp, data *new_datap )
{
 if (( *datalistpp ) == NULL ) {
     new_datap->next           = NULL;
     new_datap->prev           = NULL;
     (*datalistpp)                = new_datap;
 } else {
     new_datap->next                = (*datalistpp);
     new_datap->prev         = NULL;
     (*datalistpp)->prev        = new_datap;
     (*datalistpp)                = new_datap;
 }
}


void    calculate( data *data_list )
{
 data        *one, *two;
 double        cor, cor2;
 int                i;

 /* First, change all the data in datalist[] to be subtraction of each's mean */
 for ( one = data_list; one != NULL; one = one->next ) {
     for ( i = 0; i < one->numdata; i++ ) {
         one->datalist[i]        -= one->mean;
     }
 }


 for ( one = data_list; one != NULL; one = one->next ) {
     for ( two = one->next; two != NULL; two = two->next ) {
         cor        = 0.0;

         for ( i = 0; i < one->numdata; i++ ) {
             cor        += one->datalist[i] * two->datalist[i];
         }
         cor                /= ( one->subpart * two->subpart );
         cor2        = cor * cor;

         fprintf( stdout,
                  "%s\t%s\t%lf\t%lf\n",
                  one->element_name,
                  two->element_name,
                  cor, cor2 );
     }
 }
}

void        remove_trailing_newline( char *linebuffer )
{
 size_t        len;

 while( 1 ) {
     len                = strlen( linebuffer );
     if ( linebuffer[len-1] == '\n' ) {
         linebuffer[len-1]        = '\0';
     } else {
         break;
     }
 }
}


void
help( char *programname )
{
 fprintf( stderr,
          "%s <filename containing target>\n",
          programname );
 exit( 127 );
}

void*
o_malloc( size_t size, const char *msg, const char *func, const long line )
{
 void        *p;
 p        = malloc( size );
 if ( p == NULL ) {
     fprintf( stderr, "%s:%ld: %s", func, line, msg );
     exit( 1 );
 }
}

みての通り、ファイルを読んでくるところの処理が延々。実際の計算部分は calculate() の部分だけ。しかもこれすら標準出力への出力コードが含まれているのに、このサイズ。あぁ、Cの最大の弱点はファイルIOの記述/メモリ管理の記述が面倒くさいことだなぁ。

ちなみに。「C++の例外処理で書けばいいじゃない」という提案をいただいた。が、その手段は使えないし、趣味として使わない。例外処理は素晴らしい機能だ。ただし、あれは「コード的には正しいが、動作環境との兼ね合いで不幸に至った場合」に、どこで不幸を食らったかに依存せず「不幸を食らった事実だけ」を書けばよいから素晴らしいのだ。

しかし、上記のコード、多分端々に見えているだろうが、「おかしくなるのは私のコーディングが間違っているから」というケースに対応するためのものだ。実際に発生した現象としては:

1. tagname を引っ張ってきたのはいいが、行末の "\n" を削り忘れていた。
   
2. Correlation のOK を case insensitive にチェックしているが、実は最初のコードでは "OK" に対して入力は "OK\n" で合致しない攻撃をくらいまくっていた。
   
3. Correlation が NG の場合、計算対象からはずすためにリンクリストに登録しないのだが、そのメモリを「放置」していた(メモリリークだ)。64bit環境では特に問題は無かった(ちょっと swap IO するだけ)のだが、32bit環境に持ってきたらパンクした。
   
4. Correlation 自身はOKであっても、データが全フィールドに渡って存在するとは限らない。存在しないエントリは 0.0 として扱うことになっていたのだが、sscanf() の戻り値をチェックするときにそのことを完全に失念していて、おかしな部分で give up して exit() していた。
   

まぁ、そういうケースに対応するための「デバッガコード」であって、異常系ではなかったりする。例外処理だと、「どこでおかしいのか」も判らないし、特に 1 のケースなどは assert() にもひっかけようがない。

先にも書いたが、1データ組み合わせ当たり 6.74usec で計算できた。はーやーいー。9msec とは比べ物にならない。やっぱりコンパイラ言語は本気を出すとすごいわ。

リファクタリング-4-

リファクタリング-3-で書いた「データ列を1データ列1ファイルに分割する」プログラムをまず、作った。

remove_NaN.pl:

#! /usr/bin/perl
#
# remove_NaN.pl <dataFileName> <template of eachoutput>
# chop dataFile into each row. If each row is clear that
# "Pearson product-moment correlation coefficient" calculation
# will ends in NaN, do not output that entry.
#
# For <template of eachoutput> I mean something like "hogehoge%08d".
# output file will have row number as for %08d part (%08d can be %d
# or anything in integer format ).

# use Math::NumberCruncher;
# $NCref = Math::NumberCruncher->new();

our $at_once        = 7000;

$datafilename        = $ARGV[0];
$outtemplate        = $ARGV[1];

if ( $outtemplate eq '' ) {
  $outtemplate        = 'temp%08d.TSV';
}

open( DATAFILE, $datafilename ) or die "unable to open $datafilename as datafile\n";


# First, read first line. It's TAG-NAME list.

our @tagname;
our $line;

$line                = <DATAFILE>;
@tagname        = split /\t/, $line;
for ( $i = 0; $i < scalar( @tagname ); $i++ ) {
  $tagname[$i]        =~ s/^\s*(".*")\s*$/\1/;
}


# Okey. Now let's read each rows.
# First row is not numeral, but is time of when data was sampled. so skip first row.

our        $filecounter        = 0;
for ( $row = 0; $row < scalar( @tagname ); $row += $at_once ) {
  debug( "$row start\n" );
  # bring file pointer to top of line.
  seek( DATAFILE, 0, 0 );

  # skip tagname line.
  $line        = <DATAFILE>;

  my @output;
  my @data;
  my @sum;
  my @mean;
  my @sqsum;
  my @subpart;
  my @datas;

  for ( my $i = 0; $i < $at_once; $i++ ) {
      $output[$i]        = 1;
  }

  while (<DATAFILE>) {
      my @oneline        = split /\t/, $_;
      my $val;

      for ( my $i = 0; $i < $at_once; $i++ ) {
          $oneline[$row + $i]        =~ s/^\s*"(.*)"\s*$/\1/;
          $val        = $oneline[$row + $i];
          if ( $val eq '' ) {
              $output[$i]        = 0;
          } else {
              push @{$data[$i]}, $val;
              $sum[$i]        += $val;
          }
      }
  }
  debug( "read data\n" );

  # calculate $mean and $subpart.

  for ( my $i = 0; $i < $at_once; $i++ ) {
      if ( $output[$i] ) {
          $datas[$i]        = scalar( @{$data[$i]} );
          $mean[$i]        = $sum[$i] / $datas[$i];
          if ( $mean[$i] eq 'NaN' ) {
              $output[$i]        = 0;
              next;
          }

          for ( my $j = 0; $j < $datas[$i]; $j++ ) {
              my $val;
              $val        = $data[$i][$j] - $mean[$i];
              $sqsum[$i]        += $val * $val;
          }
          $subpart[$i]        = sqrt( $sqsum[$i] );

          if (( $subpart[$i] eq 'NaN' )||( $subpart[$i] == 0.0 )) {
              $output[$i]        = 0;
          }
      }
  }
  debug( "calc data\n" );

  # dumpout the row.
  for ( my $i = 0; $i < $at_once; $i++ ) {
      my $filename        = sprintf( $outtemplate, $filecounter );
      if ( !defined( $tagname[$row+$i] )){
          next;
      }

      open( OUTPUT, "> $filename" ) or die "unable to open $filename as output file. Template was $outtemplate.\n";
      print OUTPUT "tagname\t$tagname[$row+$i]\n";
      print OUTPUT "Correlation\t" . ( $output[$i] ? 'OK' : 'NG' ) . "\n";
      print OUTPUT "sum\t$sum[$i]\n";
      print OUTPUT "mean\t$mean[$i]\n";
      print OUTPUT "sqsum\t$sqsum[$i]\n";
      print OUTPUT "subpart\t$subpart[$i]\n";
      print OUTPUT "datas\t$datas[$i]\n";
      print OUTPUT "\n";

      foreach ( @{$data[$i]} ) {
          print OUTPUT "$_\n";
      }

      close OUTPUT;
      debug( "output $filename \n" );

      $filecounter++;
  }

ROWSTEP_CLEANUP:
  undef @output;
  undef @data;
  undef @sum;
  undef @mean;
  undef @subpart;
  undef @datas;
}

close DATAFILE;

sub debug {
  my $time        = gmtime();
  print STDERR "$time: ";
  foreach ( @_ ) {
      print STDERR $_;
  }
}

このプログラム中、唯一判らないのは「$at_once」だろう。これはようするに、読んだデータの内、$at_once 個のデータ列づつ平均値だのなんだのを計算して、各データ列ごとに別のファイルに出力する、という事だ。

実はこのプログラム、かなりメモリを消費する。データファイル自体は１行づつ読んではパースし、不要な部分は捨てているのだが、元のデータが 27000列とかだと、一日分のデータだけでも 32bit Windows システムではメモリが足りずに途中でハングする。そこで、一度に処理する行数をスクリプト中に書いておくことにした。メモリが足りないならここを調節すればよい。

ただし、可能な限り大きな値を指定しておくことをお勧めする。実験として50を選んであるが、ものすごく遅くて難儀した。27000の場合、540回同じファイルを読み直すことになり、そのたびにファイル全体を文字列処理する必要があるのだから、当たり前と言えば当たり前だが。

リファクタリング-3-

Math::NumberCruncherのコードの一部を手元に展開した状態で読んでいく…ん? んん?? んんん～～～???

どうみても、Correlation() の計算式は

ピアソン積率相関係数 のそれではないっ!!!

確かに似たような値が出てきてはいるが、明らかに違うっ!!!

ぬぅ。なんということだ。大前提からしてやり直しかっ!!!

しかも。これは致命的な問題ではないからいいようなものの…同モジュールに含まれている平方根を出す関数 SqrRoot() は

0.0 を入れると NaN を返すっ!!!

まてぃぃいいいいい!!!

どこの世界に0.0 の平方根が計算できない数値演算ライブラリがあるっ!!! ちなみに、他の値を入れると(プラスの値であれば)妥当な出力を出してくる。どうやら、0.0 の場合におかしくなるような漸近式を使っているらしい。

まぁ、0.0だろうが NaN だろうが、その結果を分母に持ってきて分数を計算すれば、NaNになるのは変わらないので、そりゃ結果だけ見れば構わないと言えば構わないのだが…。とりあえず、この程度のデバッグもできていないモジュールである以上、これは信頼するわけには行かぬ。

というわけで根こそぎ作り変えてくれよう。

!!決定!!

とりあえず、次のようなデザインにする。

1. まず、データ列を1データ列1ファイルに分割する。
2. そのデータ列には、当該データ列の平均値、平均値と各データの二乗和、さらに二乗和の平方根、の3つは最低限でも、データとして記録する。
3. できればデータの個数も記録する。
4. もちろん、ラベル名も記録する。
5. いうまでもなく、データ列の値そのものも1行1エレメントで記録する。
6. 2 で得られた値のどれかが NaN になるか、あるいは「二乗和の平方根」が 0.0 になる場合は、そのようなデータ列は当該データファイルを作らない。
7. このファイルはすべて文字列で数値を記録する。

で、このファイルを複数読み込んで、組み合わせ毎に計算するプログラムを別途作る。各データ列ごとの事前計算は全部記録されているので、本来の計算に集中できる…に違いない。

…それにしてもここまで巻き戻るとは…

リファクタリング-2-

リファクタリングとして数多くのことをする必要があるが、とりあえず最初にやるべきことは、Math::NumberCruncher のコードから、今回必要なルーチンをコピーして持ってくることだろう。Correlation()関数を呼び出す回数は O(n²) になるが、Wikipediaの「相関係数」のページによると相関係数を求める式はこのようになっている:

分母が特に「列」単位で独立して計算できる事がわかる。列単位ならば計算はO(n)で済む。このような最適化がかけられるかどうか、調べるには元のソースが無くてはいけない。
幸い、CPANには perl のソースが .pm ファイルとしてついてくる。それをコピーしてみたが、全部 perl で書かれているようだ。これならば最適かもかけられよう。

次に、x と y のラベル名を読み込むコードを変更する。旧来:

our $i = 0;
our @x_tagname = ();
while (<COMPARE_X>) {
 $_ =~ s/^\"//;
 $_ =~ s/\"\n$//;
 $x_tagname[$i] = $1;
 $i++;
}
close COMPARE_X;

となっていたのだが、

1. push()関数を使えばインデックスの $i は不要になる。


2. ラベル名に " がついたままでもいいじゃないか。どうせ x, y だけでなくデータのラベルも " がついたままなのだから。文字列比較の際に " がついたまま比較すればいい

と考えるた。結果、このように簡単になる:

our @x_tagname = ();
while (<COMPARE_X>) {
 $_ =~ s/^\s*(".*")\s*$/\1/;
 push @x_tagname, ( $_ );
}
close COMPARE_X;

同じ変更は COMPARE_Y を読み込む部分についても行う。DATAFILEの最初の1行目を読んだ部分に関してはもっと極端だ。

# Read first line of .
# This should be the list of TAG-NAME.
$line  = <DATAFILE>;
@tagname = split /\t/, $line;
for ( $i = 0; $i < scalar( @tagname ); $i++ ) {
  my $val;
  $val  = $tagname[$i];
  $val  =~ s/^\s*"//;
  $val  =~ s/"\s*$//;
  $tagname[$i] = $val;
}

というこれは

$line  = <DATAFILE>;
@tagname = split /\t/, $line;
for ( my $i = 0; $i < scalar( @tagname ); $i++ ) {
    $tagname[$i] =~ s/^\s*(".*")\s*$/\1/g;
}

とできる。

次に。今回のリファクタリングで、BestFit() を呼び出し、回帰直線の傾きと切片を求めるのをやめる。この計算も結構な高コストである上に、相関係数を求めた結果として、重要と思われたものについてのみ必要な数値だからだ。これらは、計算すべきターゲットが見つかってから計算するのでも間に合う。