#include <iostream>
#include <vector>
#include <Eigen/Dense>

int main()
{
  std::vector<std::array<float, 3> > cloud;
  cloud.push_back({ -1.841658711, 0.324050009, 1.041258931 });
  cloud.push_back({-1.449076772, 0.431593835, 0.987673104 });
  cloud.push_back({-1.04681015, 0.415599823, 1.049820662 });
  cloud.push_back({-0.646572888, 0.42549479, 1.087714672 });
  cloud.push_back({-0.25920701, 0.599574149, 0.971796691 });
  cloud.push_back({0.149769783, 0.497985959, 1.114130974 });
  cloud.push_back({-1.883069754, 0.6849733, 1.249212146 });
  cloud.push_back({-1.481496215, 0.677821219, 1.303076267 });
  cloud.push_back({-1.083532453, 0.716715276, 1.313803315 });

  Eigen::MatrixXf mat(9,3);
  for (size_t i = 0; i < cloud.size();i++) {
    mat.row(i) << cloud[i][0] , cloud[i][1] , cloud[i][2];
  }

  // センタリングしたデータを作成
  Eigen::MatrixXf centered = mat.rowwise() - mat.colwise().mean();

  std::cout << "---" << std::endl;

  // センタリングされたデータに対して、列ごとに標準偏差を求め、各列の要素を標準偏差で割る
  for (int c = 0; c < 3; c++) {
    auto col = mat.col(c).array();
    double std_dev = sqrt((col - col.mean()).square().sum() / col.size());
    centered.col(c) /= std_dev;
  }

  
  std::cout << std::endl << centered << std::endl;// 表示
  std::cout << std::endl << std::endl;

  //////////////////////
  // 分散共分散行列
  Eigen::MatrixXf cov = (centered.adjoint() * centered) / double(centered.rows());

  std::cout << cov << std::endl;// 表示
  return 0;


}

  // センタリングしたデータを作成
  Eigen::MatrixXf centered = mat.rowwise() - mat.colwise().mean();

Eigen::MatrixXf mat(2, 3);

mat.row(0) << 1, 3, 5;
mat.row(1) << 2, 4, 6;

std::cout << mat.colwise().mean() << std::endl;

auto col = mat.col(c).array();

double std_dev = sqrt(
  (col/*各要素から*/ - col.mean()/*全要素の平均を引き*/)
    .square()/*それぞれ二乗して*/
      .sum()/*全て合計して*/
        / col.size()/*要素数で割り*/
)/*平方根をとる*/;

centered.col(c) /= std_dev;

//Eigen::SelfAdjointEigenSolver<Eigen::MatrixXf> es(cov);
Eigen::EigenSolver<Eigen::MatrixXf> es(cov);
Eigen::MatrixXcf eigen_vectors = es.eigenvectors();
Eigen::VectorXcf eigen_values = es.eigenvalues();

#include <iostream>
#include <vector>
#include <array>

using vector3f = std::vector<std::array<float, 3>>;

//! @brief データの平均を求める
//! @param [in] vec データの配列
//! @return 各要素の平均
std::array<float, 3> Average(const vector3f& vec) {

  // 初期化
  std::array<float, 3> ave;
  ave[0] = 0;
  ave[1] = 0;
  ave[2] = 0;

  // 各要素平均
  for (size_t i = 0; i < vec.size(); i++) {
    ave[0] += vec[i][0];
    ave[1] += vec[i][1];
    ave[2] += vec[i][2];
  }
  ave[0] /= vec.size();
  ave[1] /= vec.size();
  ave[2] /= vec.size();

  return ave;

}

//分散
std::array<float, 3> Distributed(const std::array<float, 3>& average, const vector3f& vec) {
  int n = vec.size();
  std::array<float, 3> s2{ 0,0,0 };
  for (int j = 0; j < 3; j++) {
    for (const auto& v : vec) {
      s2[j] += std::pow(v[j] - average[j], 2);
    }
    s2[j] = 1.f / n * s2[j];
  }
  return s2;
}