/**
   SOM ( Self-Organizing Maps ; 自己組織化写像 )
**/
class SOM
{
  typedef std::vector< double > Coef; // 係数ベクトル(ニューロン)の型

  Matrix< Coef > outLayer_; // 出力層
  size_t dimension_;        // 係数の次元
  unsigned int cnt_;        // 入力回数
  unsigned int maxCnt_;     // 最大入力回数

  bool hTorus_;             // 水平方向に連結
  bool vTorus_;             // 垂直方向に連結

  /*
    勝ちニューロン(Winner)の位置を求める

    in 入力層の係数ベクトル
    row,col 求める位置を保持する変数へのポインタ
  */
  void getWinner( const Coef& in, size_t* row, size_t* col );

public:

  /*
    出力層の大きさ ( row, col ) と係数ベクトルの次元 dim、係数の最大値 max、最大学習回数 maxCnt を指定して構築
  */
  SOM( size_t row, size_t col, size_t dim, double max, unsigned int maxCnt, bool hTorus, bool vTorus )
    { init( row, col, dim, max, maxCnt, hTorus, vTorus ); }

  /*
    init : 出力層を初期化する

    出力層は乱数で初期化され、学習回数もリセットされる。

    row,col 出力層の大きさ
    dim 係数ベクトルの次元
    max 係数の最大値(乱数生成の最大値として利用する)
    maxCnt 最大学習回数
  */
  void init( size_t row, size_t col, size_t dim, double max, unsigned int maxCnt, bool hTorus, bool vTorus );

  /*
    training : 学習を行う

    in 入力する係数ベクトル
    nFunc 近傍関数

    戻り値 : 学習回数が最大学習回数以上なら false を返す
  */
  bool training( const std::vector< double >& in, const SOM_NeighborhoodFunction& nFunc );

  /*
    operator[] : 出力層の行 row の定数反復子を返す

    SOM オブジェクト som から som[r][c][i] で r 行 c 列の i 番目の要素へアクセスできる
  */
  Matrix< std::vector< double > >::const_iterator operator[]( size_t row ) const
  { return( outLayer_.row( row ) ); }

  // rows : 出力層の行数を返す
  size_t rows() const
  { return( outLayer_.rows() ); }

  // cols : 出力層の列数を返す
  size_t cols() const
  { return( outLayer_.cols() ); }

  // coefSize : 係数ベクトルのサイズを返す
  size_t coefSize() const
  { return( dimension_ ); }
};

/*
  SOM::init : SOM の初期化

  row,col : 出力層の大きさ
  dim : 係数ベクトルの次元
  max : 係数の最大値(乱数生成の最大値として利用する)
  maxCnt : 最大学習回数
*/
void SOM::init( size_t row, size_t col, size_t dim, double max, unsigned int maxCnt, bool hTorus, bool vTorus )
{
  assert( row > 0 && col > 0 && dim > 0 );

  // 出力層を乱数でリセット
  dimension_ = dim;
  outLayer_.assign( row, col );
  for ( Matrix< Coef >::iterator out = outLayer_.begin() ; out != outLayer_.end() ; ++out ) {
    out->resize( dimension_ );
    for ( Coef::iterator coef = out->begin() ; coef != out->end() ; ++coef )
      *coef = static_cast< double >( rand() ) * max / ( static_cast< double >( RAND_MAX ) + 1 );
  }

  // 学習回数のリセット
  cnt_ = 0;
  maxCnt_ = maxCnt;

  // 水平・垂直方向に連結しているか
  hTorus_ = hTorus;
  vTorus_ = vTorus;
}

/*
  SOM::getWinner : 勝ちニューロンの位置を取得する

  in に最も近いニューロンの位置を ( row, col ) に返す
*/
void SOM::getWinner( const Coef& in, size_t* row, size_t* col )
{
  assert( in.size() == dimension_ );

  double min = std::numeric_limits< double >::max();

  for ( Matrix< Coef >::size_type r = 0 ; r < outLayer_.rows() ; ++r ) {
    for ( Matrix< Coef >::size_type c = 0 ; c < outLayer_.cols() ; ++c ) {
      double d = L2Distance()( in.begin(), in.end(), outLayer_[r][c].begin() );
      if ( d < min ) {
        *row = r;
        *col = c;
        min = d;
      }
    }
  }
}

/*
  GetMinRange : 更新範囲の最小値を決定する

  winner : 勝ちニューロン
  range : 更新範囲の幅(winner からの距離)
  isTorus : 端を連結するなら true
*/
int GetMinRange( size_t winner, size_t range, bool isTorus )
{
  if ( range <= winner || isTorus )
    return( static_cast< int >( winner ) - range );

  // 範囲外で端を連結させない場合はゼロを返す
  return( 0 );
}

/*
  GetMinRange : 更新範囲の最大値を決定する

  winner : 勝ちニューロン
  range : 更新範囲の幅(winner からの距離)
  size : 出力層の大きさ
  isTorus : 端を連結するなら true

  戻り値 : 学習回数が最大学習回数以上なら false を返す
*/
int GetMaxRange( size_t winner, size_t range, size_t size, bool isTorus )
{
  if ( winner + range < size || isTorus )
    return( winner + range );

  // 範囲外で端を連結させない場合は最大値を返す
  return( size - 1 );
}

/*
  toSizeT : int 型の更新範囲 i を size_t 型の SOM の位置に変換する

  size : 出力層の大きさ
*/
size_t toSizeT( int i, size_t size )
{
  while ( i < 0 )
    i += size;

  return( i % size );
}

/*
  SOM::training : 学習を行う

  in : 入力する係数ベクトル
  nFunc : 近傍関数

  戻り値 : 学習回数が最大学習回数以上なら false を返す
*/
bool SOM::training( const Coef& in, const SOM_NeighborhoodFunction& nFunc )
{
  assert( &in != 0 && &nFunc != 0 );

  if ( cnt_ >= maxCnt_ ) return( false );

  size_t wr, wc;

  getWinner( in, &wr, &wc ); // 勝ちニューロン(BMU)の位置を取得する

  // 近傍関数の値が有効な範囲をチェックする
  size_t range = 0;
  while ( nFunc( 0, 0, range, 0, cnt_, maxCnt_ ) > 1E-9 )
    ++range;

  // 更新範囲を決定
  int minR = GetMinRange( wr, range, vTorus_ );
  int minC = GetMinRange( wc, range, hTorus_ );
  int maxR = GetMaxRange( wr, range, outLayer_.rows(), vTorus_ );
  int maxC = GetMaxRange( wc, range, outLayer_.cols(), hTorus_ );

  // 出力層の更新
  for ( int mr = minR ; mr <= maxR ; ++mr ) {
    for ( int mc = minC ; mc <= maxC ; ++mc ) {
      Coef& out = outLayer_[toSizeT( mr, outLayer_.rows() )][toSizeT( mc, outLayer_.cols() )];
      for ( size_t i = 0 ; i < out.size() ; ++i ) {
        double d = nFunc( wr, wc, mr, mc, cnt_, maxCnt_ ) * ( in[i] - out[i] );
        out[i] += d;
      }
    }
  }
  ++cnt_;

  return( true );
}

/**
   SOM 専用近傍関数基底クラス
**/
struct SOM_NeighborhoodFunction
{
  /*
    operator() : 近傍関数の値を返す

    wr,wc 勝ちニューロン(BMU)の位置
    mr,mc 更新するニューロンの位置
    cnt 学習回数
    maxCnt 最大学習回数
  */
  virtual double operator()( double wr, double wc, double mr, double mc, unsigned int cnt, unsigned int maxCnt ) const = 0;

  // 仮想デストラクタ(何もしない)
  virtual ~SOM_NeighborhoodFunction() {}
};

/**
   @brief L2 ノルムを使った SOM 専用近傍関数
**/
class SOM_L2_NeighborhoodFunction : public SOM_NeighborhoodFunction
{
  double c_; // スケール調整のための定数

public:

  // スケール調整のための定数 c を指定して構築
  SOM_L2_NeighborhoodFunction( double c )
    : c_( c ) { assert( c_ > 0 ); }

  // 近傍関数の値を返す
  virtual double operator()( double wr, double wc, double mr, double mc, unsigned int cnt, unsigned int maxCnt ) const;
};

/*
  SOM_L2_NeighborhoodFunction::operator() : 近傍関数の値を返す

  wr,wc : 勝ちニューロン(BMU)の位置
  mr,mc : 更新するニューロンの位置
  cnt : 学習回数
  maxCnt : 最大学習回数
*/
double SOM_L2_NeighborhoodFunction::operator()( double wr, double wc, double mr, double mc, unsigned int cnt, unsigned int maxCnt ) const
{
  if ( cnt == maxCnt )
    return( 0 );

  return( c_ * std::exp( -( ( wr - mr ) * ( wr - mr ) + ( wc - mc ) * ( wc - mc ) ) / std::pow( 1.0 - static_cast< double >( cnt ) / maxCnt, 2.0 ) ) );
}