template< class X >
using item_type = std::map< X, unsigned int >;   // アイテムの型
template< class X >
using list_type = std::vector< item_type< X > >; // アイテムリストの型
using stat_type = double;                        // 数量の型
using size_type = size_t;                        // 添字の型

/*
  GetCategory : カテゴリ毎のリストを作成する

  [ s, e ) : カテゴリデータ
  category : 生成するカテゴリリストへのポインタ
*/
template< class In, class X >
void GetCategory( In s, In e, item_type< X, unsigned int >* category )
{
  for ( ; s != e ; ++s ) {
    auto i = category->find( *s );
    if ( i == category->end() )
      category->insert( std::make_pair( *s, 1 ) );
    else
      ++( i->second ); // カテゴリに属する要素数をカウント
  }
}

/*
  GetCategoryList : カテゴリ毎のリストを作成する

  x : アイテムごとのカテゴリデータ ( 列:アイテム 行:要素 )
  category : 生成するカテゴリリストへのポインタ
*/
template< class X >
void GetCategoryList( const Matrix< X >& x, list_type< X, unsigned int >* category )
{
  using size_type = typename Matrix< X >::size_type;

  for ( size_type c = 0 ; c < x.cols() ; ++c ) {
    category->push_back( item_type< X, unsigned int >() );
    auto& map = category->back();
    GetCategory( x.column( c ), x.column( c ).end(), &map );
  }
}

/*
  GetDupMatched : 指定した二つのアイテム・カテゴリを両方持つ要素の数をカウントする

  x : アイテムごとのカテゴリデータ ( 列:アイテム 行:要素 )
  j, u : アイテムの番号 ( x の列番号 )
  k, v : カテゴリ

  戻り値 : 計算したカウント
*/
template< class X >
stat_type GetDupMatched( const Matrix< X >& x,
                         typename Matrix< X >::size_type j, const X& k,
                         typename Matrix< X >::size_type u, const X& v )
{
  using size_type = typename Matrix< X >::size_type;

  // 同一アイテム内で異なるカテゴリが重複する事はない
  if ( j == u && k != v ) return( 0 );

  stat_type total = 0;
  for ( size_type r = 0 ; r < x.rows() ; ++r )
    if ( x[r][j] == k && x[r][u] == v )
      total += 1;

  return( total );
}

/*
  GetMatchedYSum : 指定したアイテム・カテゴリを持つ要素の従属変数の和を求める

  x : アイテムごとのカテゴリデータ ( 列:アイテム 行:要素 )
  ys : 従属変数のリストの先頭
  j : アイテムの番号 ( x の列番号 )
  k : カテゴリ

  戻り値 : 計算した和
*/
template< class X, class In >
stat_type GetMatchedYSum( const Matrix< X >& x, In ys,
                          typename Matrix< X >::size_type j, const X& k )
{
  using size_type = typename Matrix< X >::size_type;

  stat_type total = 0;
  for ( size_type r = 0 ; r < x.rows() ; ++r ) {
    if ( x[r][j] == k ) total += *ys;
    ++ys;
  }

  return( total );
}

/*
  CreateLES : 数量化 I 類用の正規方程式の係数計算

  x : アイテムごとのカテゴリデータ ( 列:アイテム 行:要素 )
  ys : 従属変数のリストの先頭
  category : カテゴリのリスト
  coef : 生成する係数行列へのポインタ
  a : 生成する解ベクトルへのポインタ
*/
template< class X, class In >
void CreateLES( const Matrix< X >& x, In ys, const list_type< X, unsigned int >& category,
                SquareMatrix< stat_type >* coef, std::vector< stat_type >* a )
{
  using SM = SquareMatrix< stat_type >;
  using List = list_type< X, unsigned int >;

  SM buff; // 線形独立の形にする前の係数行列

  // カテゴリの総数を求める
  typename List::size_type total = 0;
  for ( auto i = category.begin() ; i != category.end() ; ++i )
    total += i->size();

  // 配列・行列のリサイズ
  buff.assign( total, total );
  a->resize( total - category.size() + 1 );
  coef->assign( total - category.size() + 1, total - category.size() + 1 );

  // 係数の計算
  SM::size_type startR = 0;
  SM::size_type c = 0;
  for ( typename List::size_type j = 0 ; j < category.size() ; ++j ) {
    for ( auto k = category[j].begin() ; k != category[j].end() ; ++k ) {
      SM::size_type r = startR;
      typename List::size_type u = j;
      for ( auto v = k ; v != category[j].end() ; ++v ) {
        buff[r][c] = buff[c][r] = GetDupMatched( x, j, k->first, u, v->first );
        ++r;
      }
      while ( ++u < category.size() ) {
        for ( auto v = category[u].begin() ; v != category[u].end() ; ++v ) {
          buff[r][c] = buff[c][r] = GetDupMatched( x, j, k->first, u, v->first );
          ++r;
        }
      }
      ++c;
      ++startR;
    }
  }

  // 解ベクトルの計算
  auto ai = a->begin();
  for ( typename List::size_type j = 0 ; j < category.size() ; ++j ) {
    for ( auto k = category[j].begin() ; k != category[j].end() ; ++k ) {
      if ( j > 0 && k == category[j].begin() ) continue; // 各アイテムの最初のカテゴリは除外(最初のアイテムを除く)
      *ai++ = GetMatchedYSum( x, ys, j, k->first );
    }
  }

  // 線形独立の形になるように行列を削除する
  std::set< SM::size_type > ignoreRC; // 削除する行列番号
  SM::size_type i = 0;
  for ( typename List::size_type j = 0 ; j + 1 < category.size() ; ++j ) {
    i += category[j].size();
    ignoreRC.insert( i );
  }
  SM::size_type oR = 0; // 出力行
  for ( SM::size_type r = 0 ; r < total ; ++r ) {
    if ( ignoreRC.find( r ) != ignoreRC.end() ) continue;
    SM::size_type oC = 0; // 出力列
    for ( SM::size_type c = 0 ; c < total ; ++c ) {
      if ( ignoreRC.find( c ) != ignoreRC.end() ) continue;
      (*coef)[oR][oC] = buff[r][c];
      ++oC;
    }
    ++oR;
  }
}

/*
  NormalizeCoef : 係数の標準化

  x : アイテムごとのカテゴリデータ ( 列:アイテム 行:要素 )
  category : カテゴリのリスト
  a :  標準化する係数
*/
template< class X >
void NormalizeCoef( const Matrix< X >& x, const list_type< X, unsigned int >& category,
                    list_type< X, stat_type >* a )
{
  using size_type = typename list_type< X, unsigned int >::size_type;

  for ( size_type j = 0 ; j < category.size() ; ++j ) {
    // 減算する値の計算 Δ = Σ( a_jk * n_jk ) / N
    stat_type delta = 0;
    auto& aj = (*a)[j];
    for ( auto k = category[j].begin() ; k != category[j].end() ; ++k )
      delta += aj[k->first] * k->second;
    delta /= x.rows();
    // 標準化
    for ( auto k = category[j].begin() ; k != category[j].end() ; ++k )
      aj[k->first] -= delta;
  }
}

/*
  CalcSjj : x の分散値を求める

  x : アイテムごとのカテゴリデータ ( 列:アイテム 行:要素 )
  j : アイテム番号
  category : カテゴリのリスト
  aj : 対象アイテムの係数

  戻り値 : 求めた分散値
*/
template< class X >
stat_type CalcSjj( const Matrix< X >& x, typename Matrix< X >::size_type j,
                   const list_type< X, unsigned int >& category,
                   const item_type< X, stat_type >& aj )
{
  stat_type t1 = 0;
  stat_type t2 = 0;
  for ( auto k = category[j].begin() ; k != category[j].end() ; ++k ) {
    auto ajk = aj.find( k->first );
    t1 += std::pow( ajk->second, 2 ) * k->second; // Σ( a_jk * a_jk * n_jk )
    t2 += ajk->second * k->second;                // Σ( a_jk * n_jk )
  }

  return( ( t1 - std::pow( t2, 2 ) / x.rows() ) / x.rows() );
}

/*
  CalcSyj : x と y の共分散を求める

  x : アイテムごとのカテゴリデータ ( 列:アイテム 行:要素 )
  j : アイテム番号
  category : カテゴリのリスト
  ys : y のリストの先頭
  my : y の平均値
  aj : 対象アイテムの係数

  戻り値 : 求めた共分散値
*/
template< class X, class In >
stat_type CalcSyj( const Matrix< X >& x, typename Matrix< X >::size_type j,
                   const list_type< X, unsigned int >& category, In ys, stat_type my,
                   const item_type< X, stat_type >& aj )
{
  stat_type t1 = 0;
  stat_type t2 = 0;
  for ( auto k = category[j].begin() ; k != category[j].end() ; ++k ) {
    auto ajk = aj.find( k->first );
    t1 += ajk->second * GetMatchedYSum( x, ys, j, k->first ); // Σ( a_jk * Σy_iδ_ijk )
    t2 += ajk->second * k->second;                            // Σ( a_jk * n_jk )
  }

  return( ( t1 - my * t2 ) / x.rows() );
}

/*
  CalcSju : x 同士の共分散を求める

  x : アイテムごとのカテゴリデータ ( 列:アイテム 行:要素 )
  j, u : アイテム番号
  category : カテゴリのリスト
  aj, au : 対象アイテムの係数

  戻り値 : 求めた共分散値
*/
template< class X >
stat_type CalcSju( const Matrix< X >& x,
                   typename Matrix< X >::size_type j, typename Matrix< X >::size_type u,
                   const list_type< X, unsigned int >& category,
                   const item_type< X, stat_type >& aj, const item_type< X, stat_type >& au )
{
  const auto& cj = category[j];
  const auto& cu = category[u];

  stat_type t1 = 0;
  for ( auto cjk = cj.begin() ; cjk != cj.end() ; ++cjk ) {
    for ( auto cuv = cu.begin() ; cuv != cu.end() ; ++cuv ) {
      auto ajk = aj.find( cjk->first );
      auto auv = au.find( cuv->first );
      // Σ( a_jk * a_uv * Σ( δ_ijk * δ_iuv ) )
      t1 += ajk->second * auv->second * GetDupMatched( x, j, cjk->first, u, cuv->first );
    }
  }

  stat_type t2 = 0;
  for ( auto cjk = cj.begin() ; cjk != cj.end() ; ++cjk ) {
    auto ajk = aj.find( cjk->first );
    t2 += ajk->second * cjk->second; // Σ( a_jk * n_jk )
  }

  stat_type t3 = 0;
  for ( auto cuv = cu.begin() ; cuv != cu.end() ; ++cuv ) {
    auto auv = au.find( cuv->first );
    t3 += auv->second * cuv->second; // Σ( a_uv * n_uv )
  }

  return( ( t1 - t2 * t3 / x.rows() ) / x.rows() );
}

/*
  CreateCovariance : 共分散行列の計算

  x : アイテムごとのカテゴリデータ ( 列:アイテム 行:要素 )
  category : カテゴリのリスト
  ys, ye : y の先頭と終端の次の位置
  my : y の平均値
  a : 係数
  s : 共分散行列を保持する変数のポインタ
*/
template< class X, class In >
void CreateCovariance( const Matrix< X >& x, const list_type< X, unsigned int >& category,
                       In ys, In ye, stat_type my, const list_type< X, stat_type >& a,
                       SquareMatrix< stat_type >* s )
{
  using size_type = SquareMatrix< stat_type >::size_type;

  s->assign( x.cols() + 1 );

  // syy の計算
  (*s)[0][0] = Statistics::SampleVariance< stat_type >( ys, ye );

  // sjj の計算
  for ( size_type j = 1 ; j < s->rows() ; ++j )
    (*s)[j][j] = CalcSjj( x, j - 1, category, a[j - 1] );

  // syj の計算
  for ( size_type j = 1 ; j < s->rows() ; ++j )
    (*s)[j][0] = (*s)[0][j] = CalcSyj( x, j - 1, category, ys, my, a[j - 1] );

  // sju の計算
  for ( size_type j = 1 ; j < s->rows() ; ++j )
    for ( size_type u = j + 1 ; u < s->rows() ; ++u )
      (*s)[j][u] = (*s)[u][j] = CalcSju( x, j - 1, u - 1, category, a[j - 1], a[u - 1] );
}

/*
  数量化 I 類
*/
template< class X >
class Quantification1
{
  using x_item_type = item_type< X, unsigned int >; // 独立変数用のアイテムの型
  using x_list_type = list_type< X, unsigned int >; // 独立変数用のアイテムリストの型
  using a_item_type = item_type< X, stat_type >;    // カテゴリごとの係数の型
  using a_list_type = list_type< X, stat_type >;    // 全アイテムの係数リストの型

  a_list_type a_;                  // 係数
  x_list_type category_;           // カテゴリリスト
  stat_type my_;                   // y の平均値
  SquareMatrix< stat_type > s_;    // 共分散行列
  SquareMatrix< stat_type > invS_; // 共分散行列の逆行列

public:

  /*
    独立変数と従属変数を指定して構築

    独立変数ベクトル(X)は、一行に一つのベクトルの要素が、列ごとに各ベクトルが並んだ行列で表される。
    ベクトルの数(標本数)は従属変数(Y)の要素数で決まる。そのため、X行列の行数とYの要素数は一致
    させる必要がある。

    独立ベクトルの各要素はカテゴリ(分類)データとして扱われ、内部で同一データごとに集計される。
    データの型は文字列・数値データのいずれも利用できるが、浮動小数点数などは少しでも値が異なれば
    違うデータとして扱われるため注意すること。

    x : 独立変数(X)を保持する行列
    ys, ye : 従属変数(Y)のデータ列先頭と末尾の次
    solver : 連立方程式を求めるソルバー
  */
  template< class In, class Solver >
  Quantification1( const Matrix< X >& x, In ys, In ye, Solver solver );

  /*
    アイテム数を返す
  */
  size_type itemSize() const
  { return( category_.size() ); }

  /*
    指定したアイテムのカテゴリ数を返す
  */
  size_type categorySize( size_type item ) const
  { return( category_.at( item ).size() ); }

  /*
    回帰係数 ( 数量 ) を返す

    アイテム数や、各アイテムのカテゴリ数より大きな添字を指定した場合は std::out_of_range 例外を投げる。

    item : アイテムの番号
    category : カテゴリの番号
  */
  stat_type a( size_type item, const X& category ) const;

  /*
    カテゴリ名を返す

    アイテム数や、各アイテムのカテゴリ数より大きな添字を指定した場合は std::out_of_range 例外を投げる。

    item : アイテムの番号
    category : カテゴリの番号
  */
  X category( size_type item, size_type category ) const;

  /*
    従属変数の予測値を返す

    xs, xe : 独立変数の先頭と末尾の次の位置
  */
  template< class In >
  stat_type expY( In xs, In xe ) const;

  /*
    重相関係数を返す

  */
  stat_type mcc() const
  { return( std::sqrt( 1.0 - 1.0 / ( s_[0][0] * invS_[0][0] ) ) ); }

  /*
    指定したアイテム番号と y の間の偏相関係数を返す
  */
  stat_type partialCC( size_type item ) const;
};

/*
  Quantification1 コンストラクタ

  x : アイテムごとのカテゴリデータ ( 列:アイテム 行:要素 )
  ys, ye : y の先頭と終端の次の位置
  solver : 連立方程式を解くためのソルバ
*/
template< class X >
template< class In, class Solver >
Quantification1< X >::Quantification1( const Matrix< X >& x, In ys, In ye, Solver solver )
{
  ErrLib::CheckLinearModel( x, ys, ye );

  // カテゴリリストの作成
  GetCategoryList( x, &category_ );

  // y の平均値を求める
  my_ = Statistics::Mean< stat_type >( ys, ye );

  // 正規方程式の作成
  SquareMatrix< stat_type > coef;
  std::vector< stat_type > aBuff;
  CreateLES( x, ys, category_, &coef, &aBuff );

  // 係数の計算
  solver( coef, aBuff.begin() );

  // 足りない係数(ゼロとおいた係数)を加える
  auto j = category_.begin();
  auto a = aBuff.begin();
  a_.push_back( a_item_type() );
  for ( auto k = j->begin() ; k != j->end() ; ++k )
    ( a_.back() )[k->first] = *a++;
  for ( ++j ; j != category_.end() ; ++j ) {
    auto k = j->begin();
    a_.push_back( a_item_type() );
    ( a_.back() )[k->first] = 0;
    for ( ++k ; k != j->end() ; ++k )
      ( a_.back() )[k->first] = *a++;
  }

  // 係数の標準化
  NormalizeCoef( x, category_, &a_ );

  // 共分散行列の計算
  CreateCovariance( x, category_, ys, ye, my_, a_, &s_ );
  MathLib::InverseMatrix( s_, &invS_, solver );
}

/*
  Quantification1< X >::a : 指定したアイテム番号・カテゴリに対する係数 ( 数量 ) を返す

  item : アイテム番号
  category : カテゴリ
*/
template< class X >
typename Quantification1< X >::stat_type
Quantification1< X >::a( size_type item, const X& category ) const
{
  const std::string FUNC_NAME = "Quantification1< X >::a";

  // アイテム番号が範囲外なら out_of_range 例外を投げる
  if ( item >= category_.size() )
    throw std::out_of_range( FUNC_NAME + " : item index is out of range" );

  // カテゴリが見つからなかったら out_of_range 例外を投げる
  auto i = a_[item].find( category );
  if ( i == a_[item].end() )
    throw std::out_of_range( FUNC_NAME + " : category not found" );

  return( i->second );
}

/*
  Quantification1< X >::category : 指定したアイテム・カテゴリ番号のカテゴリの値を返す

  item : アイテム番号
  category : カテゴリ番号

  戻り値 : カテゴリの値
*/
template< class X >
X Quantification1< X >::category( size_type item, size_type category ) const
{
  const std::string FUNC_NAME = "Quantification1< X >::category";

  // 範囲外なら out_of_range 例外を投げる
  if ( item >= category_.size() )
    throw std::out_of_range( FUNC_NAME + " : item index is out of range" );
  if ( category >= category_[item].size() )
    throw std::out_of_range( FUNC_NAME + " : category index is out of range" );

  const auto& categoryList = category_[item];
  auto ans = categoryList.begin();
  std::advance( ans, category );

  return( ans->first );
}

/*
  Quantification1< X >::expY : 予測値を返す

  xs xe : 独立変数の開始と終了の次の位置

  戻り値 : 予測値
*/
template< class X >
template< class In >
typename Quantification1< X >::stat_type
Quantification1< X >::expY( In xs, In xe ) const
{
  const std::string FUNC_NAME = "Quantification1< X >::expY";

  ErrLib::CheckArrayLength( xs, xe, itemSize() );

  stat_type ans = my_;
  for ( size_type item = 0 ; xs != xe ; ++xs, ++item ) {
    auto ajk = a_[item].find( *xs );
    if ( ajk == a_[item].end() )
      throw std::out_of_range( FUNC_NAME + " : category not found" );
    ans += ajk->second;
  }

  return( ans );
}

/*
  Quantification1< X >::partialCC : 指定したアイテム番号と y の間の偏相関係数を返す

  item : アイテム番号
*/
template< class X >
typename Quantification1< X >::stat_type
Quantification1< X >::partialCC( size_type item ) const
{
  const std::string FUNC_NAME = "Quantification1< X >::partialCC";

  // アイテム番号が範囲外なら out_of_range 例外を投げる
  if ( item >= category_.size() )
    throw std::out_of_range( FUNC_NAME + " : item index is out of range" );

  return( -invS_[0][item + 1] / std::sqrt( invS_[0][0] * invS_[item + 1][item + 1] ) );
}

/*
  CalcNijk : あるカテゴリに属する要素の中で、指定したアイテム・カテゴリの独立変数を持つものをカウントする

  x : 独立変数
  ys : 従属変数の開始位置
  j : 対象のアイテム番号
  cjk : 対象のカテゴリ
  y : 要素の対象カテゴリ

  戻り値 : カウントした結果
*/
template< class X, class In >
stat_type CalcNijk( const Matrix< X >& x, In ys,
                    typename Matrix< X >::size_type j, const X& cjk,
                    typename std::iterator_traits< In >::value_type y )
{
  using size_type = typename Matrix< X >::size_type;

  stat_type total = 0;
  for ( size_type r = 0 ; r < x.rows() ; ++r ) {
    if ( *ys++ != y ) continue;
    if ( x[r][j] == cjk )
      total += 1;
  }

  return( total );
}

/*
  CalcFjk_uv : f( jk, uv ) の計算

  x : 独立変数
  j, u : 対象のアイテム番号
  k, v : 対象のカテゴリ反復子

  戻り値 : f( jk, uv ) の値
*/
template< class X >
stat_type CalcFjk_uv( const Matrix< X >& x,
                      typename list_type< X, unsigned int >::size_type j, typename item_type< X, unsigned int >::const_iterator k,
                      typename list_type< X, unsigned int >::size_type u, typename item_type< X, unsigned int >::const_iterator v )
{
  stat_type t1 = GetDupMatched( x, j, k->first, u, v->first ); // Σ( δ_ijk * δ_iuv )
  stat_type t2 = k->second * v->second / x.rows();             // n_jk * n_uv / N

  return( t1 - t2 );
}

/*
  CalcGjk_uv : g( jk, uv ) の計算

  x : 独立変数
  ys : 従属変数の開始位置
  j, u : 対象のアイテム番号
  k, v : 対象のカテゴリ反復子
  y : 要素のカテゴリリスト

  戻り値 : g( jk, uv ) の値
*/
template< class X, class In >
stat_type CalcGjk_uv( const Matrix< X >& x, In ys,
                      typename list_type< X, unsigned int >::size_type j, typename item_type< X, unsigned int >::const_iterator k,
                      typename list_type< X, unsigned int >::size_type u, typename item_type< X, unsigned int >::const_iterator v,
                      const item_type< typename std::iterator_traits< In >::value_type, unsigned int >& y )
{
  stat_type t1 = 0;
  for ( auto i = y.begin() ; i != y.end() ; ++i ) {
    stat_type n_ijk = CalcNijk( x, ys, j, k->first, i->first );
    stat_type n_iuv = CalcNijk( x, ys, u, v->first, i->first );
    t1 += n_ijk * n_iuv / i->second; // Σ( n_ijk * n_iuv / n_i )
  }
  stat_type t2 = k->second * v->second / x.rows(); // n_jk * n_uv / N

  return( t1 - t2 );
}

/*
  CreateMatrix : 行列 F, G の計算

  x : 独立変数
  ys : 従属変数の開始位置
  xCategory : カテゴリリスト(独立変数用)
  yCategory : カテゴリリスト(従属変数用)
  f : 行列 F へのポインタ
  g : 行列 G へのポインタ
*/
template< class X, class In >
void CreateMatrix( const Matrix< X >& x, In ys,
                   const list_type< X, unsigned int >& xCategory,
                   const item_type< typename std::iterator_traits< In >::value_type, unsigned int >& yCategory,
                   SymmetricMatrix< stat_type >* f, SymmetricMatrix< stat_type >* g )
{
  using SM = SymmetricMatrix< stat_type >;
  using List = list_type< X, unsigned int >;

  SM fBuff, gBuff; // 線形独立の形にする前の係数行列

  // カテゴリの総数を求める
  typename List::size_type total = 0;
  for ( auto i = xCategory.begin() ; i != xCategory.end() ; ++i )
    total += i->size();

  // 行列のリサイズ
  fBuff.assign( total, total );
  gBuff.assign( total, total );
  f->assign( total - xCategory.size(), total - xCategory.size() );
  g->assign( total - xCategory.size(), total - xCategory.size() );

  // 行列の要素を求める
  typename List::size_type j = 0;
  auto k = xCategory[j].begin();
  for ( SM::size_type r = 0 ; r < fBuff.rows() ; ++r ) {
    typename List::size_type u = j;
    auto v = k;
    for ( SM::size_type c = r ; c < fBuff.cols() ; ++c ) {
      fBuff[r][c] = CalcFjk_uv( x, j, k, u, v );
      gBuff[r][c] = CalcGjk_uv( x, ys, j, k, u, v, yCategory );
      if ( ++v == xCategory[u].end() ) {
        ++u;
        v = xCategory[u].begin();
      }
    }
    if ( ++k == xCategory[j].end() ) {
      ++j;
      k = xCategory[j].begin();
    }
  }

  // 線形独立の形になるように行列を削除する
  std::set< SM::size_type > ignoreRC; // 削除する行列番号
  SM::size_type i = 0;
  for ( typename List::size_type j = 0 ; j < xCategory.size() ; ++j ) {
    ignoreRC.insert( i );
    i += xCategory[j].size();
  }
  SM::size_type oR = 0; // 出力行
  for ( SM::size_type r = 0 ; r < total ; ++r ) {
    if ( ignoreRC.find( r ) != ignoreRC.end() ) continue;
    SM::size_type oC = oR; // 出力列
    for ( SM::size_type c = r ; c < total ; ++c ) {
      if ( ignoreRC.find( c ) != ignoreRC.end() ) continue;
      (*f)[oR][oC] = fBuff[r][c];
      (*g)[oR][oC] = gBuff[r][c];
      ++oC;
    }
    ++oR;
  }
}

/*
  GeneralEigenvalueProblem : 一般固有値問題 ( G - λF )a = 0 を解く

  f : 行列 F
  g : 行列 G
  r : 固有値を保持する行列へのポインタ (対角成分が固有値になる)
  q : 固有ベクトルを保持する行列へのポインタ
  e : 収束判定のためのしきい値
  maxCnt : 収束しなかった場合の最大処理回数
  solver : 連立方程式の計算オブジェクト ( 逆行列計算用 )
*/
template< class Solver >
void GeneralEigenvalueProblem( const SymmetricMatrix< stat_type >& f, const SymmetricMatrix< stat_type >& g,
                               SquareMatrix< stat_type >* r, SquareMatrix< stat_type >* q,
                               stat_type e, unsigned int maxCnt, Solver solver )
{
  // F を固有値分解
  try {
    EigenLib::Householder_DoubleShiftQR( f, r, q, e, maxCnt );
  } catch( ExceptionExcessLimit& ee ) {
    SymmetricMatrix< stat_type > eigen( f );
    EigenLib::Jacobi( &eigen, q, e, maxCnt );
    *r = eigen;
  }

  // P = Q√R の計算
  SquareMatrix< stat_type > p;
  p.assign( r->rows() );
  for ( SquareMatrix< stat_type >::size_type i = 0 ; i < r->rows() ; ++i )
    for ( SquareMatrix< stat_type >::size_type j = 0 ; j < q->rows() ; ++j )
      p[j][i] = (*q)[j][i] * std::sqrt( (*r)[i][i] );

  // P の転置行列
  SquareMatrix< stat_type > pt( p );
  pt.transpose();

  // P, Pt の逆行列
  SquareMatrix< stat_type > invP, invPt;
  MathLib::InverseMatrix( p, &invP, solver );
  MathLib::InverseMatrix( pt, &invPt, solver );

  // P^-1GPt^-1 の計算
  SquareMatrix< stat_type > buff;
  MatrixLib::mult( invP, g, &buff );
  MatrixLib::mult( buff, invPt, &buff );

  SymmetricMatrix< stat_type > pmGPtm;
  pmGPtm.assign( buff.size() );
  for ( SymmetricMatrix< stat_type >::size_type i = 0 ; i < buff.rows() ; ++i )
    for ( SymmetricMatrix< stat_type >::size_type j = i ; j < buff.cols() ; ++j )
      pmGPtm[i][j] = buff[i][j];

  // P^-1GPt^-1 を固有値分解
  try {
    EigenLib::Householder_DoubleShiftQR( pmGPtm, r, q, e, maxCnt );
  } catch( ExceptionExcessLimit& ee ) {
    EigenLib::Jacobi( &pmGPtm, q, e, maxCnt );
    *r = pmGPtm;
  }
  for ( SquareMatrix< stat_type >::size_type j = 0 ; j < q->cols() ; ++j )
    MatrixLib::multMV( invPt, q->column( j ), q->column( j ) );
}

/*
  CalcBi : 従属変数のカテゴリごとに Y の平均を求める

  x : 独立変数
  ys : 従属変数の開始位置
  yCategory : カテゴリリスト(従属変数用)
  a : 求めた係数
  bi : Y の平均を保持する配列へのポインタ
*/
template< class X, class In >
void CalcBi( const Matrix< X >& x,
             In ys, const item_type< typename std::iterator_traits< In >::value_type, unsigned int >& yCategory,
             const list_type< X, stat_type >& a,
             item_type< typename std::iterator_traits< In >::value_type, stat_type >* bi )
{
  // bi の初期化
  for ( auto i = yCategory.begin() ; i != yCategory.end() ; ++i )
    bi->insert( std::make_pair( i->first, 0 ) );

  // Y の和をカテゴリごとに求める
  for ( typename Matrix< X >::size_type c = 0 ; c < x.cols() ; ++c ) {
    const auto& ac = a[c];
    In y = ys;
    for ( typename Matrix< X >::size_type r = 0 ; r < x.rows() ; ++r ) {
      auto i = ac.find( x[r][c] );
      (*bi)[*y++] += i->second;
    }
  }

  // カテゴリごとの個数で割って平均とする
  for ( auto i = bi->begin() ; i != bi->end() ; ++i ) {
    auto y = yCategory.find( i->first );
    i->second /= y->second;
  }
}

/*
  CalcSyy : y の分散値を求める ( 数量化 II 類専用 )

  x : 独立変数
  xCategory : カテゴリリスト(独立変数用)
  ys : 従属変数の開始位置
  yCategory : カテゴリリスト(従属変数用)
  a : 求めた係数値
  bi : 従属変数のカテゴリごとの y の平均
*/
template< class X, class In >
stat_type CalcSyy( const Matrix< X >& x, const list_type< X, unsigned int >& xCategory,
                   In ys, const item_type< typename std::iterator_traits< In >::value_type, unsigned int >& yCategory,
                   const list_type< X, stat_type >& a,
                   const item_type< typename std::iterator_traits< In >::value_type, stat_type >& bi )
{
  stat_type t1 = 0;
  stat_type t2 = 0;
  for ( auto i = bi.begin() ; i != bi.end() ; ++i ) {
    auto j = yCategory.find( i->first );
    t1 += std::pow( i->second, 2 ) * j->second; // n_i * b_i * b_i
    t2 += i->second * j->second;                // n_i * b_i
  }

  return( ( t1 - std::pow( t2, 2 ) / x.rows() ) / x.rows() );
}

/*
  CalcSyj : x と y の共分散を求める

  x : アイテムごとのカテゴリデータ ( 列:アイテム 行:要素 )
  j : アイテム番号
  xCategory : カテゴリのリスト ( 独立変数用 )
  ys : y のリストの先頭
  yCategory : カテゴリのリスト ( 従属変数用 )
  aj : 対象アイテムの係数
  bi : 従属変数のカテゴリごとの y の平均
*/
template< class X, class In >
stat_type CalcSyj( const Matrix< X >& x, typename Matrix< X >::size_type j,
                   const list_type< X, unsigned int >& xCategory,
                   In ys, const item_type< typename std::iterator_traits< In >::value_type, unsigned int >& yCategory,
                   const item_type< X, stat_type >& aj,
                   const item_type< typename std::iterator_traits< In >::value_type, stat_type >& bi )
{
  stat_type t1 = 0;
  stat_type t2 = 0;
  for ( auto k = xCategory[j].begin() ; k != xCategory[j].end() ; ++k ) {
    stat_type sum = 0;
    for ( auto i = bi.begin() ; i != bi.end() ; ++i )
      sum += i->second * CalcNijk( x, ys, j, k->first, i->first );
    auto ajk = aj.find( k->first );
    t1 += sum * ajk->second;       // ΣkΣi( a_jk * b_i * n_ijk )
    t2 += ajk->second * k->second; // Σk( a_jk * n_jk )
  }
  stat_type t3 = 0;
  for ( auto i = yCategory.begin() ; i != yCategory.end() ; ++i ) {
    t3 += ( bi.find( i->first ) )->second * i->second; // Σi( n_i * b_i )
  }

  return( ( t1 - t2 * t3 / x.rows() ) / x.rows() );
}

/*
  CreateCovariance : 共分散行列の計算

  x : アイテムごとのカテゴリデータ ( 列:アイテム 行:要素 )
  xCategory : カテゴリのリスト ( 独立変数用 )
  ys : y のリストの先頭
  yCategory : カテゴリのリスト ( 従属変数用 ) 
  a : 求めた係数
  bi : 従属変数のカテゴリごとの y の平均
  s : 共分散行列を保持する変数のポインタ
*/
template< class X, class In >
void CreateCovariance( const Matrix< X >& x, const list_type< X, unsigned int >& xCategory,
                       In ys, const item_type< typename std::iterator_traits< In >::value_type, unsigned int >& yCategory,
                       const list_type< X, stat_type >& a, const item_type< typename std::iterator_traits< In >::value_type, stat_type >& bi,
                       SquareMatrix< stat_type >* s )
{
  using size_type = SquareMatrix< stat_type >::size_type;

  s->assign( x.cols() + 1 );

  // syy の計算
  (*s)[0][0] = CalcSyy( x, xCategory, ys, yCategory, a, bi );

  // sjj の計算
  for ( size_type j = 1 ; j < s->rows() ; ++j )
    (*s)[j][j] = CalcSjj( x, j - 1, xCategory, a[j - 1] );

  // syj の計算
  for ( size_type j = 1 ; j < s->rows() ; ++j )
    (*s)[j][0] = (*s)[0][j] = CalcSyj( x, j - 1, xCategory, ys, yCategory, a[j - 1], bi );

  // sju の計算
  for ( size_type j = 1 ; j < s->rows() ; ++j )
    for ( size_type u = j + 1 ; u < s->rows() ; ++u )
      (*s)[j][u] = (*s)[u][j] = CalcSju( x, j - 1, u - 1, xCategory, a[j - 1], a[u - 1] );
}

/*
  数量化 II 類
*/
template< class X, class Y >
class Quantification2
{
  using x_item_type = item_type< X, unsigned int >; // 独立変数用のアイテムの型
  using x_list_type = list_type< X, unsigned int >; // 独立変数用のアイテムリストの型
  using y_type = item_type< Y, unsigned int >;      // 従属変数用のカテゴリリストの型
  using bi_type = item_type< Y, stat_type >;        // カテゴリごとの y 値の型
  using a_item_type = item_type< X, stat_type >;    // カテゴリごとの係数の型
  using a_list_type = list_type< X, stat_type >;    // 全アイテムの係数リストの型

  a_list_type a_;                  // 係数
  x_list_type xCategory_;          // x のアイテムリスト
  y_type yCategory_;               // y のカテゴリリスト
  bi_type bi_;                     // 予測値(yに対する数量)
  SquareMatrix< stat_type > s_;    // 共分散行列
  SquareMatrix< stat_type > invS_; // 共分散行列の逆行列

public:

  /*
    独立変数と従属変数を指定して構築

    独立変数ベクトル(X)は、一行に一つのベクトルの要素が、列ごとに各ベクトルが並んだ行列で表される。
    ベクトルの数(標本数)は従属変数(Y)の要素数で決まる。そのため、X行列の行数とYの要素数は一致
    させる必要がある。

    独立ベクトルの各要素はカテゴリ(分類)データとして扱われ、内部で同一データごとに集計される。
    データの型は文字列・数値データのいずれも利用できるが、浮動小数点数などは少しでも値が異なれば
    違うデータとして扱われるため注意すること。

    x : 独立変数(X)を保持する行列
    ys, ye : 従属変数(Y)のデータ列先頭と末尾の次
    solver : 連立方程式を求めるソルバー
  */
  template< class In, class Solver >
  Quantification2( const Matrix< X >& x, In ys, In ye, Solver solver );

  /*
    アイテム数を返す
  */
  size_type itemSize() const
  { return( xCategory_.size() ); }

  /*
    指定したアイテムのカテゴリ数を返す

    アイテム番号が範囲外だった場合、例外 std::out_of_range を投げる。
  */
  size_type categorySize( size_type item ) const
  { return( xCategory_.at( item ).size() ); }

  /*
    回帰係数(数量)を返す

    アイテム数や、各アイテムのカテゴリ数より大きな添字を指定した場合は std::out_of_range 例外を投げる。

    item : アイテムの番号
    category : カテゴリ
  */
  stat_type a( size_type item, const X& category ) const;

  /*
    指定したアイテム・カテゴリ番号のカテゴリの値を返す

    item : アイテム番号
    category : カテゴリ番号
  */
  X category( size_type item, size_type category ) const;

  /*
    指定した独立変数から予測値(数量)を求める

    s, e : 独立変数の開始と末尾の次の位置
  */
  template< class In >
  stat_type expY( In s, In e ) const;

  /*
    指定した独立変数から判別を行う

    s, e : 独立変数の開始と末尾の次の位置
  */
  template< class In >
  Y discriminant( In s, In e ) const;

  /*
    重相関係数を返す
  */
  stat_type mcc() const
  { return( std::sqrt( 1.0 - 1.0 / ( s_[0][0] * invS_[0][0] ) ) ); }

  /*
    指定したアイテム番号と y の間の偏相関係数を返す
  */
  stat_type partialCC( size_type item ) const;
};

/*
  Quantification2 コンストラクタ

  x : アイテムごとのカテゴリデータ ( 列:アイテム 行:要素 )
  ys, ye : y の先頭と終端の次の位置
  solver : 連立方程式を解くためのソルバ
*/
template< class X, class Y >
template< class In, class Solver >
Quantification2< X, Y >::Quantification2( const Matrix< X >& x, In ys, In ye, Solver solver )
{
  const stat_type THRESHOLD = 1E-6;
  const unsigned int MAX_CNT = 10000;

  const std::string FUNC_NAME = "Quantification2 constructor";

  ErrLib::CheckLinearModel( x, ys, ye );

  GetCategoryList( x, &xCategory_ );
  GetCategory( ys, ye, &yCategory_ );

  SymmetricMatrix< stat_type > f, g;
  CreateMatrix( x, ys, xCategory_, yCategory_, &f, &g );

  SquareMatrix< stat_type > r, q;
  GeneralEigenvalueProblem( f, g, &r, &q, THRESHOLD, MAX_CNT, solver );

  // 最大固有値の固有ベクトルを抽出する
  if ( r.size() == 0 )
    throw std::runtime_error( "No eigenvalue found" );
  stat_type maxEigenValue = r[0][0];
  SquareMatrix< stat_type >::size_type maxEigenVector = 0;
  for ( SquareMatrix< stat_type >::size_type i = 1 ; i < r.size() ; ++i ) {
    if ( maxEigenValue < r[i][i] ) {
      maxEigenValue = r[i][i];
      maxEigenVector = i;
    }
  }
  SquareMatrix< stat_type >::const_iterator c = q.column( maxEigenVector );
  for ( auto i = xCategory_.begin() ; i != xCategory_.end() ; ++i ) {
    auto j = i->begin();
    a_.push_back( a_item_type() );
    ( a_.back() )[j->first] = 0;
    for ( ++j ; j != i->end() ; ++j )
      ( a_.back() )[j->first] = *c++;
  }

  // 係数の標準化
  NormalizeCoef( x, xCategory_, &a_ );

  CalcBi( x, ys, yCategory_, a_, &bi_ );

  CreateCovariance( x, xCategory_, ys, yCategory_, a_, bi_, &s_ );
  MathLib::InverseMatrix( s_, &invS_, solver );
}

/*
  Quantification2< X, Y >::a : 指定したアイテム番号・カテゴリに対する係数 ( 数量 ) を返す

  item : アイテム番号
  category : カテゴリ
*/
template< class X, class Y >
typename Quantification2< X, Y >::stat_type
Quantification2< X, Y >::a( size_type item, const X& category ) const
{
  const std::string FUNC_NAME = "Quantification2< X, Y >::a";

  // アイテム番号が範囲外なら out_of_range 例外を投げる
  if ( item >= xCategory_.size() )
    throw std::out_of_range( FUNC_NAME + " : item index is out of range" );

  // カテゴリが見つからなかったら out_of_range 例外を投げる
  auto i = a_[item].find( category );
  if ( i == a_[item].end() )
    throw std::out_of_range( FUNC_NAME + " : category not found" );

  return( i->second );
}

/*
  Quantification2< X, Y >::category : 指定したアイテム・カテゴリ番号のカテゴリの値を返す

  item : アイテム番号
  category : カテゴリ番号

  戻り値 : カテゴリの値
*/
template< class X, class Y >
X Quantification2< X, Y >::category( size_type item, size_type category ) const
{
  const std::string FUNC_NAME = "Quantification2< X, Y >::category";

  // 範囲外なら out_of_range 例外を投げる
  if ( item >= xCategory_.size() )
    throw std::out_of_range( FUNC_NAME + " : item index is out of range" );
  if ( category >= xCategory_[item].size() )
    throw std::out_of_range( FUNC_NAME + " : category index is out of range" );

  const auto& categoryList = xCategory_[item];
  auto ans = categoryList.begin();
  std::advance( ans, category );

  return( ans->first );
}

/*
  Quantification2< X, Y >::expY : 独立変数から y の予測値を返す

  s, e : 独立変数の先頭と末尾の次の位置
*/
template< class X, class Y >
template< class In >
typename Quantification2< X, Y >::stat_type
Quantification2< X, Y >::expY( In s, In e ) const
{
  const std::string FUNC_NAME = "Quantification2< X, Y >::expY";

  ErrLib::CheckArrayLength( s, e, a_.size() );

  stat_type y = 0;
  for ( typename a_list_type::size_type j = 0 ; j < a_.size() ; ++j ) {
    const auto& aj = a_[j];
    const auto& ajk = aj.find( *s++ );
    if ( ajk == aj.end() )
      throw std::out_of_range( FUNC_NAME + " : category not found" );
    y += ajk->second;
  }

  return( y );
}

/*
  Quantification2< X, Y >::discriminant : 独立変数から y のカテゴリを判別する

  s, e : 独立変数の先頭と末尾の次の位置
*/
template< class X, class Y >
template< class In >
Y Quantification2< X, Y >::discriminant( In s, In e ) const
{
  const std::string FUNC_NAME = "Quantification2< X, Y >::discriminant";

  // y の予測値に最も近い bi を探す
  stat_type ey = expY( s, e );
  auto minIndex = bi_.begin();
  if ( minIndex == bi_.end() )
    return( ProbDist::nan() );
  stat_type minDiff = std::abs( minIndex->second - ey );
  auto i = minIndex;
  while ( ++i != bi_.end() ) {
    stat_type diff = std::abs( i->second - ey );
    if ( minDiff > diff ) {
      minDiff = diff;
      minIndex = i;
    }
  }

  return( minIndex->first );
}

/*
  Quantification2< X, Y >::partialCC : 指定したアイテム番号と y の間の偏相関係数を返す

  item : アイテム番号
*/
template< class X, class Y >
typename Quantification2< X, Y >::stat_type
Quantification2< X, Y >::partialCC( size_type item ) const
{
  const std::string FUNC_NAME = "Quantification2< X, Y >::partialCC";

  // アイテム番号が範囲外なら out_of_range 例外を投げる
  if ( item >= xCategory_.size() )
    throw std::out_of_range( FUNC_NAME + " : item index is out of range" );

  return( -invS_[0][item + 1] / std::sqrt( invS_[0][0] * invS_[item + 1][item + 1] ) );
}

/*
  CalcFi : 要素 i が反応するカテゴリの数 fi を求める

  data : 対象データ
  fi : 求めた fi を保持する配列へのポインタ
*/
void CalcFi( const Matrix< bool >& data, std::vector< stat_type >* fi )
{
  fi->clear();

  for ( Matrix< bool >::size_type i = 0 ; i < data.rows() ; ++i ) {
    stat_type total = 0;
    for ( Matrix< bool >::size_type j = 0 ; j < data.cols() ; ++j )
      if ( data[i][j] ) total += 1;
    fi->push_back( total );
  }
}

/*
  CalcGj : カテゴリ j に反応した要素の数 gj を求める

  data : 対象データ
  gj : 求めた gj を保持する配列へのポインタ
*/
void CalcGj( const Matrix< bool >& data, std::vector< stat_type >* gj )
{
  gj->clear();

  for ( Matrix< bool >::size_type j = 0 ; j < data.cols() ; ++j ) {
    stat_type total = 0;
    for ( Matrix< bool >::size_type i = 0 ; i < data.rows() ; ++i )
      if ( data[i][j] ) total += 1;
    gj->push_back( total );
  }
}

/*
  CreateMatrixJ : yj に対する固有値問題の対称行列を計算する

  data : 対象データ
  fi : 要素 i が反応するカテゴリの数
  gj : カテゴリ j に反応した要素の数
  mat : 求めた行列を保持する変数へのポインタ
*/
void CreateMatrixJ( const Matrix< bool >& data,
                    const std::vector< stat_type >& fi, const std::vector< stat_type >& gj,
                    SymmetricMatrix< stat_type >* mat )
{
  mat->assign( data.cols() );
  for ( Matrix< bool >::size_type j = 0 ; j < data.cols() ; ++j ) {
    for ( Matrix< bool >::size_type v = j ; v < data.cols() ; ++v ) {
      (*mat)[j][v] = 0;
      for ( Matrix< bool >::size_type i = 0 ; i < data.rows() ; ++i ) {
        if ( data[i][j] && data[i][v] )
          (*mat)[j][v] += 1.0 / ( fi[i] * std::sqrt( gj[j] ) * std::sqrt( gj[v] ) );
      }
    }
  }
}

/*
  CreateMatrixI : xi に対する固有値問題の対称行列を計算する

  data : 対象データ
  fi : 要素 i が反応するカテゴリの数
  gj : カテゴリ j に反応した要素の数
  mat : 求めた行列を保持する変数へのポインタ
*/
void CreateMatrixI( const Matrix< bool >& data,
                    const std::vector< stat_type >& fi, const std::vector< stat_type >& gj,
                    SymmetricMatrix< stat_type >* mat )
{
  mat->assign( data.rows() );
  for ( Matrix< bool >::size_type i = 0 ; i < data.rows() ; ++i ) {
    for ( Matrix< bool >::size_type u = i ; u < data.rows() ; ++u ) {
      (*mat)[i][u] = 0;
      for ( Matrix< bool >::size_type j = 0 ; j < data.cols() ; ++j ) {
        if ( data[i][j] && data[u][j] )
          (*mat)[i][u] += 1.0 / ( gj[j] * std::sqrt( fi[i] ) * std::sqrt( fi[u] ) );
      }
    }
  }
}

/*
  数量化 III 類
*/
class Quantification3
{
  std::vector< std::vector< stat_type > > xi_; // 各要素の数量
  std::vector< std::vector< stat_type > > yj_; // 各カテゴリの数量
  std::vector< stat_type > lambda_;            // 固有値

public:

  /*
    各要素のカテゴリに対する反応を指定して構築

    引数 data は、一行に一つの要素のカテゴリに対する反応が、一列に一つのカテゴリの要素ごとの反応が並んだ行列で表される。

    data : 各要素のカテゴリ毎の反応を保持する行列
  */
  Quantification3( const Matrix< bool >& data );

  /*
    要素数を返す
  */
  size_type dataSize() const;

  /*
    カテゴリ数を返す
  */
  size_type categorySize() const;

  /*
    固有値の数を返す
  */
  size_type lambdaSize() const
  { return( lambda_.size() ); }

  /*
    指定した要素の数量を返す

    固有値数・要素数より大きな添字を指定した場合は std::out_of_range 例外を投げる。

    k : 固有値の順位 ( 降順 )
    i : 要素の番号
  */
  stat_type x( size_type k, size_type i ) const
  { return( xi_.at( k ).at( i ) ); }

  /*
    指定したカテゴリの数量を返す

    固有値数・カテゴリ数より大きな添字を指定した場合は std::out_of_range 例外を投げる。

    k : 固有値の順位 ( 降順 )
    j : カテゴリの番号
  */
  stat_type y( size_type k, size_type j ) const
  { return( yj_.at( k ).at( j ) ); }

  /*
    指定した順位の固有値を返す

    固有値の数より大きな添字を指定した場合は std::out_of_range 例外を投げる。

    k : 固有値の順位 ( 降順 )
  */
  stat_type lambda( size_type k ) const
  { return( lambda_.at( k ) ); }
};

/*
  Quantification3 コンストラクタ : データ data を指定して構築
*/
Quantification3::Quantification3( const Matrix< bool >& data )
{
  const stat_type THRESHOLD = 1E-6;   // 固有値計算用しきい値
  const unsigned int MAX_CNT = 10000; // 固有値計算用最大反復回数

  // fi, gj の計算
  vector< stat_type > fi, gj;

  CalcFi( data, &fi );
  CalcGj( data, &gj );

  // 固有値分解を行う行列の計算
  SymmetricMatrix< stat_type > mat;
  if ( data.cols() < data.rows() )
    CreateMatrixJ( data, fi, gj, &mat );
  else
    CreateMatrixI( data, fi, gj, &mat );

  // 固有値分解
  SquareMatrix< stat_type > r, q;
  try {
    EigenLib::Householder_DoubleShiftQR( mat, &r, &q, THRESHOLD, MAX_CNT );
  } catch( ExceptionExcessLimit& ee ) {
    SymmetricMatrix< stat_type > eigen( mat );
    EigenLib::Jacobi( &eigen, &q, THRESHOLD, MAX_CNT );
    r = eigen;
  }

  // 固有値の取得
  std::multimap< stat_type, SquareMatrix< stat_type >::size_type, std::greater< stat_type > > l; // 固有値用バッファ
  for ( SquareMatrix< stat_type >::size_type i = 0 ; i < r.rows() ; ++i )
    l.insert( std::make_pair( r[i][i], i ) );
  l.erase( l.begin() ); // 最大値(=1)は無効解なので除外

  // xi, yj の取得
  if ( data.cols() < data.rows() ) {
    for ( auto lp = l.begin() ; lp != l.end() ; ++lp ) {
      lambda_.push_back( lp->first );
      yj_.push_back( vector< stat_type >( q.column( lp->second ), q.column( lp->second ).end() ) );
      xi_.push_back( vector< stat_type >( fi.size(), 0 ) );
      auto& y = yj_.back();
      auto& x = xi_.back();
      for ( vector< stat_type >::size_type j = 0 ; j < gj.size() ; ++j )
        y[j] /= std::sqrt( gj[j] );
      for ( vector< stat_type >::size_type i = 0 ; i < fi.size() ; ++i ) {
        for ( vector< stat_type >::size_type j = 0 ; j < gj.size() ; ++j )
          if ( data[i][j] )
            x[i] += y[j] / fi[i];
        x[i] /= lp->first;
      }
    }
  } else {
    for ( auto lp = l.begin() ; lp != l.end() ; ++lp ) {
      lambda_.push_back( lp->first );
      xi_.push_back( vector< stat_type >( q.column( lp->second ), q.column( lp->second ).end() ) );
      yj_.push_back( vector< stat_type >( gj.size(), 0 ) );
      auto& x = xi_.back();
      auto& y = yj_.back();
      for ( vector< stat_type >::size_type i = 0 ; i < fi.size() ; ++i )
        x[i] /= std::sqrt( fi[i] );
      for ( vector< stat_type >::size_type j = 0 ; j < gj.size() ; ++j )
        for ( vector< stat_type >::size_type i = 0 ; i < fi.size() ; ++i ) {
          if ( data[i][j] )
            y[j] += x[i] / gj[j];
        y[j] /= lp->first;
      }
    }
  }
}

/*
  Quantification3::dataSize : 要素数を返す
*/
Quantification3::size_type Quantification3::dataSize() const
{
  if ( xi_.size() == 0 ) return( 0 );

  return( xi_[0].size() );
}

/*
  Quantification3::categorySize : カテゴリ数を返す
*/
Quantification3::size_type Quantification3::categorySize() const
{
  if ( yj_.size() == 0 ) return( 0 );

  return( yj_[0].size() );
}

/*
  CreateH : 親近性からなる行列 E から H を作成する

  e : 親近性からなる行列 E
  h : 作成した行列 H を保持する変数へのポインタ
*/
void CreateH( const SquareMatrix< stat_type >& e, SymmetricMatrix< stat_type >* h )
{
  h->assign( e.size() );

  // 非対角成分の計算
  for ( SquareMatrix< stat_type >::size_type r = 0 ; r < e.size() ; ++r ) {
    for ( SquareMatrix< stat_type >::size_type c = r + 1 ; c < e.size() ; ++c ) {
      ( *h )[r][c] = e[r][c] + e[c][r];
    }
  }

  // 対角成分の計算
  for ( SquareMatrix< stat_type >::size_type r = 0 ; r < e.size() ; ++r ) {
    stat_type total = 0;
    for ( SquareMatrix< stat_type >::size_type c = 0 ; c < e.size() ; ++c ) {
      if ( c == r ) continue;
      total += ( *h )[r][c];
    }
    ( *h )[r][r] = -total;
  }
}

/*
  数量化 IV 類
*/
class Quantification4
{
  std::vector< std::vector< stat_type > > xi_; // 各要素の数量
  std::vector< stat_type > lambda_;            // 固有値

public:

  /*
    親近性を指定して構築
  */
  Quantification4( const SquareMatrix< stat_type >& e );

  /*
    要素数を返す
  */
  size_type dataSize() const;

  /*
    固有値の数を返す
  */
  size_type lambdaSize() const
  { return( lambda_.size() ); }

  /*
    指定した数量を返す

    固有値数・要素数より大きな添字を指定した場合は std::out_of_range 例外を投げる。

    k : 固有値の順位 ( 降順 )
    i : 要素の番号
  */
  stat_type x( size_type k, size_type i ) const
  { return( xi_.at( k ).at( i ) ); }

  /*
    指定した順位の固有値を返す

    固有値の数より大きな添字を指定した場合は std::out_of_range 例外を投げる。

    k 固有値の順位 ( 降順 )
  */
  stat_type lambda( size_type k ) const
  { return( lambda_.at( k ) ); }
};

/*
  Quantification4 コンストラクタ : 親近性からなる行列 E を指定して構築

  e : 親近性からなる行列 E
*/
Quantification4::Quantification4( const SquareMatrix< stat_type >& e )
{
  const stat_type THRESHOLD = 1E-6;   // 固有値計算用しきい値
  const unsigned int MAX_CNT = 10000; // 固有値計算用最大反復回数

  // 行列 H の作成
  SymmetricMatrix< stat_type > h;
  CreateH( e, &h );

  // 固有値分解
  SquareMatrix< stat_type > r, q;
  try {
    EigenLib::Householder_DoubleShiftQR( h, &r, &q, THRESHOLD, MAX_CNT );
  } catch( ExceptionExcessLimit& ee ) {
    SymmetricMatrix< stat_type > eigen( h );
    EigenLib::Jacobi( &eigen, &q, THRESHOLD, MAX_CNT );
    r = eigen;
  }

  // 固有値の取得
  multimap< stat_type, SquareMatrix< stat_type >::size_type, std::greater< stat_type > > l; // 固有値用バッファ
  for ( SquareMatrix< stat_type >::size_type i = 0 ; i < r.rows() ; ++i )
    l.insert( std::make_pair( r[i][i], i ) );

  // 数量の取得
  for ( auto lp = l.begin() ; lp != l.end() ; ++lp ) {
    lambda_.push_back( lp->first );
    xi_.push_back( vector< stat_type >( q.column( lp->second ), q.column( lp->second ).end() ) );
  }
}

/*
  Quantification4::dataSize : 要素数を返す
*/
Quantification4::size_type Quantification4::dataSize() const
{
  if ( xi_.size() == 0 ) return( 0 );

  return( xi_[0].size() );
}