// pca.C Version 7.44   Copyleft (c) 杉本典夫
// Time-stamp: <2025-10-31 17:33:49 tac>
//
// 指定した項目について主成分分析を行う。
//
// 使用ライブラリ:libdans.so, libfun.so
// コンパイル法 :g++ -O -s -W pca pca.C
//             -I/usr/local/include -L/usr/local/lib -ldans -lfun

#include <dans.H>               // dans用ヘッダファイル(fun.H,mat.H含む)

using namespace dans;           // dans名前空間使用(名前空間matvec含む)

const vecstr USAGE =            // プログラム使用法表示用
{
  "pca  Version 7.44            Copyleft(Vo\\oV) 2025 杉本典夫",
  "書式:pca  [<オプション> …]  [<ファイル> …]",
  "機能:相関、共分散、積和行列に基づく主成分分析をしたろやないか。",
  "<オプション>",
  "  -c=<累積寄与率>:累積寄与率基準値を指定するんや。(0≦累積寄与率≦1)",
  "  -dnsname[=<出力ファイル名>]:DANSデータを出力するファイル名を指定します。",
  "    オプション省略時:DANSデータ出力せず  -dnsname:-dnsname=pca_dns.d",
  "    出力項目:主成分スコア",
  "  -f=<ファイル>:指定したファイルからオプションを読み込むんや。",
  "  -h[=<コマンド名>]:オンラインマニュアルを表示するんや。(-h:-h=man)",
  "  -m={r|s|v|x}:計算開始行列の種類を指定するんや。",
  "    r:相関行列(オプション省略時)  s:積和行列  v:共分散行列  x:単純積和行列",
  "  -ngpname[=<出力ファイル名>]:Ngraph用データの出力ファイル名を指定します。",
  "    オプション省略時:Ngraph用データ出力せず  -ngpname:-ngpname=pca_ngp.dat",
  "    出力項目:主成分負荷量",
  "  -q=<主成分数>:主成分抽出時の最大数を指定するんや。(0<主成分数≦変数の数)",
  "  -res=[-:[{def|all|pc|pcc|pcl|stat}][,…]:解析結果の出力項目を指定します。",
  "    -:-を付けた項目を出力しない(-res=-:全ての項目を出力しない)",
  "    def:オプション省略時と同じ項目を出力する  all:全ての項目を出力する",
  "    pc:主成分  pcc:主成分スコア計算用係数  pcl:主成分負荷量",
  "    stat:基礎統計量",
  "  -std[={-|+}]:主成分スコアを標準化するかどうか指定するんや。(-std:std=+)"
};

const string Title = "=== 主成分分析(PCA:principal component analysis) ==="; // 表題
const string DNSNAMEdef = "pca_dns.d"; // 省略時DANSファイル名
const string NGPNAMEdef = "pca_ngp.dat"; // 省略時Ngraphファイル名
constexpr double Cdef = 0.8;    // 省略時累積寄与率

class option                    // オプションクラス
{
private:
  string v_path;                // プログラムのパス
  string v_pname;               // プログラム名
  double v_c;                   // 累積寄与率基準値
  string v_dnsname;             // DANSファイル名
  mattype v_m;                  // 計算開始行列の種類
  string v_ngpname;             // Ngraphファイル名
  int v_q;                      // 主成分抽出用最大数
  lstr reslst;                  // 解析結果出力リスト
  mapsb resmap;                 // 解析結果出力の有無マップ
  bool v_std;                   // 主成分スコアの標準化フラッグ
public:
  option(const char* c_pname = nullptr); // コンストラクタ
  ~option()                              // デストラクタ
  {}
  const string& path() const    // プログラムのパスを返す
  { return v_path; }
  const string& pname() const   // プログラム名を返す
  { return v_pname; }
  double c() const              // 累積寄与率基準値を返す
  { return v_c; }
  const string& dnsname() const // DANSデータファイル名を返す
  { return v_dnsname; }
  mattype m() const             // 計算開始行列の種類を返す
  { return v_m; }
  const string& ngpname() const // Ngraphデータファイル名を返す
  { return v_ngpname; }
  int q() const                 // 主成分抽出用最大数を返す
  { return v_q; }
  bool res(const string& code); // 分析結果出力フラッグを返す
  bool std() const              // 主成分スコアの標準化フラッグを返す
  { return v_std; }
  void setdef()                 // 省略時オプション指定を設定する
  { getoptdef(*this); }
  void set(const string& optstr, char devc = OPTdev);
  // オプション指定によってオプションを設定する
  int check(size_t ni);         // オプションのチェックとセット
  int set_res();                // -resオプションのセット
};

option::option(const char* c_pname) : v_path(), v_pname(), v_c(Cdef), v_dnsname(), v_m(matr), v_ngpname(), v_q(0), reslst(), resmap(), v_std(true)
// コンストラクタ
{
  if (c_pname == nullptr) return; // 初期設定せず
  v_pname = getname(c_pname);     // パス名なしのプログラム名取得
  v_path = c_pname;
  strsize_t idx = v_path.find(v_pname);
  if (idx != string::npos) v_path.erase(idx); // プログラム名削除
}

bool option::res(const string& code) // 分析結果出力フラッグを返す
{
  mapsb_citr resitr = resmap.find(code);
  if (resitr == resmap.end()) return false; // code無し！
  return resmap[code];                      // true/falseを返す
}

void option::set(const string& optstr, char devc)
// オプション指定によってオプションを設定する
{
  string opt[2];                // オプションのopt[0]:項目,opt[1]:内容
  if (getopt(opt, optstr, devc) < 1) return; // オプション指定なし
  int err = 0;
  if (opt[0] == "c")            // 累積寄与率基準値の指定
  {
    v_c = opt[1].empty() ? Cdef : strtodf(opt[1]);
    if (v_c < 0.0 || 1.0 < v_c) err = badopt(pname(), optstr, USAGE, 0.0, 1.0);
  }
  else if (opt[0] == "dnsname") // DANSファイル名指定
  {
    if (opt[1].empty()) v_dnsname = DNSNAMEdef;
    else v_dnsname = opt[1];
    if (v_dnsname.size() < 3) err = badopt(pname(), optstr, USAGE);
  }
  else if (opt[0] == "f")       // オプション指定ファイル
  { err = fgetopt(*this, opt[1], true); }
  else if (opt[0] == "h")       // オンラインヘルプ表示
  { err = help(pname(), USAGE, opt[1]); }
  else if (opt[0] == "ngpname")     // Ngraphファイル名指定
  {
    if (opt[1].empty()) v_ngpname = NGPNAMEdef;
    else v_ngpname = opt[1];
    if (v_ngpname.size() < 3) err = badopt(pname(), optstr, USAGE);
  }
  else if (opt[0] == "m")       // 計算開始行列指定
  {
    if (opt[1].empty() || opt[1] == "r") v_m = matr;
    else if (opt[1] == "s") v_m = mats;
    else if (opt[1] == "v") v_m = matv;
    else if (opt[1] == "x") v_m = matetc;
    else err = badopt(pname(), optstr, USAGE);
  }
  else if (opt[0] == "q")       // 最大主成分数指定
  {
    char* ptr = nullptr;
    v_q = strtol(opt[1].c_str(), &ptr, 0);
    if (v_q < 1) err = badopt(pname(), optstr, USAGE, 1.0, (double)INT_MAX);
  }
  else if (opt[0] == "res")     // 分析結果の出力指定
  {
    if (strsplit(reslst, opt[1], CHRdelm) < 1) err = badopt(pname(), optstr, USAGE);
  }
  else if (opt[0] == "std")     // 主成分スコアの標準化指定
  {
    if (opt[1].empty() || opt[1] == "+") v_std = true;
    else if (opt[1] == "-") v_std = false;
    else err = badopt(pname(), optstr, USAGE);
  }
  else                          // 間違ったオプション指定
  { err = nonopt(pname(), optstr, USAGE); }
  if (err > 0) throw help_except("ヘルプ表示");
  if (err < 0) throw option_error(string("オプション指定違反→") + optstr);
}

int option::check(size_t ni)    // オプションのチェックとセット
{
  if (ni < 2) return disperr(pname(), "項目が2個以上必要だよ！");
  if (v_q > (int)ni) return disperr(pname(), "-qの主成分最大数が項目数より大きいよ！");
  if (set_res() < 0) return -3; // -resオプションのエラー
  return 0;                     // エラーなし
}

int option::set_res()           // -resオプションのセット
{
  resmap["res"] = resmap["def"] = resmap["pc"] = resmap["pcc"] = resmap["pcl"] = resmap["stat"] = true;
  if (reslst.size() > 0)
  {
    resmap["def"] = resmap["pc"] = resmap["pcc"] = resmap["pcl"] = resmap["stat"] = false;
    for (const string& resstr : reslst)
    {
      mapsb_citr resitr = resmap.find(resstr);
      if (resstr == "-")        // 全て出力せず
      { resmap["res"] = false; }
      else if (resstr[0] == '-' && resstr.size() > 1) // 特定の項目を出力せず
      {
        resitr = resmap.find(resstr.substr(1));
        if (resitr == resmap.end())
        { return disperr(pname(), "-res=" + resstr + "オプション指定が不適切！"); }
        else
        { resmap[resstr.substr(1)] = false; }
      }
      else if (resstr == "def" || resstr == "all")
      {
        resmap["def"] = resmap["pc"] = resmap["pcc"] = resmap["pcl"] = resmap["stat"] = true;
      }
      else if (resitr == resmap.end())
      { return disperr(pname(), "-res=" + resstr + "オプション指定が不適切！"); }
      else
      { resmap[resstr] = true; }
    }
  }
  return 0;                     // エラーなし
}

class varclass                 // 変数クラス
{
private:
  mat<ldbl> v_smat;             // 積和行列
  vecldbl v_mevec;              // 平均値ベクトル
  vecldbl v_vrvec;              // 分散ベクトル
  vecldbl v_evalvec;            // 固有値ベクトル
  mat<ldbl> v_amat;             // 固有ベクトルの行列
  int v_ni;                     // 項目数
  int v_n;                      // 例数
  int v_q;                      // 主成分数
  int v_nrep;                   // 反復計算回数
  restype v_resvar;             // 計算結果の種類
  restype cal_pca(option& opt); // 主成分を抽出する
public:
  varclass(size_t ni) : v_smat(ni + 1), v_mevec(ni + 1, 0.0L), v_vrvec(ni + 1, 0.0L), v_evalvec(ni + 1, 0.0L), v_amat(ni + 1), v_ni(ni), v_n(0), v_q(0), v_nrep(0), v_resvar(resnon)
  // コンストラクタ
  {}
  ~varclass()                   // デストラクタ
  {}
  const mat<ldbl>& smat() const // 積和行列を返す
  { return v_smat; }
  const vecldbl& mevec() const  // 平均値ベクトルを返す
  { return v_mevec; }
  const vecldbl& vrvec() const  // 分散ベクトルを返す
  { return v_vrvec; }
  const vecldbl& evalvec() const // 固有値ベクトルを返す
  { return v_evalvec; }
  const mat<ldbl>& amat() const // 固有ベクトルの行列を返す
  { return v_amat; }
  int ni() const                // 項目数を返す
  { return v_ni; }
  int n() const                 // 例数を返す
  { return v_n; }
  int q() const                 // 主成分数を返す
  { return v_q; }
  ldbl eval(int idx) const      // 固有値を返す
  { return (0 <= idx && idx <= ni()) ? v_evalvec[idx] : 0.0L; }
  ldbl sxx(int idx) const       // xの分散を返す
  { return (0 <= idx && idx <= ni()) ? v_smat(idx, idx) : 0.0L; }
  int nrep() const              // 反復計算回数を返す
  { return v_nrep; }
  restype resvar() const        // 計算結果の種類を返す
  { return v_resvar; }
  void cal(option& opt, ldata& datalst); // データを集計して各種統計量を求める
};

void varclass::cal(option& opt, ldata& datalst)
// データを集計して各種統計量を求める
{
  int dim = ni() + 1;
  vecdbl fdatavec(dim, 0.0);
  fdatavec[0] = 1.0;
  for (auto& cdatai : datalst)
  {
    bool blank = false;
    for (int j = 1; blank == false && j < dim; ++j)
    {
      fdatavec[j] = cdatai.fdata(j - 1);
      if (fdatavec[j] == BLANKf) blank = true;
    }
    if (blank)
    { cdatai.sdata(0, BLANKs); } // 後の計算のために最初の項目をブランクにする
    else
    {
      for (int i = 0; i < dim; ++i)
      {
        double x = fdatavec[i];
        for (int j = i; j < dim; ++j) v_smat.add(i, j, (ldbl)(x * fdatavec[j]));
      }
    }
  }
  for (int i = 0; i < ni(); ++i)
  {
    for (int j = i + 1; j < dim; ++j) v_smat.assign(j, i, v_smat(i, j));
  }
  v_n = (int)(v_smat(0, 0));
  mat<ldbl> tmat(v_smat);       // 一時行列
  tmat.sweep(0);                // 積和行列にする
  for (int i = 0; i < dim; ++i)
  {
    v_mevec[i] = tmat(0, i);
    v_vrvec[i] = (v_n > 1) ? tmat(i, i) / (ldbl)(v_n - 1) : 0.0L;
  }
  if (opt.m() != matetc) v_smat.sweep(0); // 積和行列にする
  if (v_n < 2) v_resvar = reserr;
  else v_resvar = cal_pca(opt);
}

restype varclass::cal_pca(option& opt)
// 主成分を抽出する
{
  ldbl w = (n() > 1) ? 1.0L / (ldbl)(n() - 1) : 0.0L;
  mat<ldbl> tmat(v_smat.sub(1, 1, ni(), ni()));
  switch (opt.m())
  {
  case matr : smat2rmat(tmat);
    break;
  case matv : tmat *= w;
  }
  v_evalvec[0] = 0.0;          // 計算開始行列のトレース(寄与率の分母)
  for (int i = 0; i < ni(); ++i) v_evalvec[0] += tmat(i, i);
  v_q = 0;                      // 主成分数
  ldbl seval = 0.0L;            // 累積寄与率
  vec<ldbl> evecv(ni());        // 固有ベクトル
  int qmax = (opt.q() > 0) ? opt.q() : ni();
  restype respca = resnon;
  do
  {
    ldbl eval = eigenud(tmat, evecv); // 固有値抽出
    if (sgn(eval) <= 0)
    { respca = resok; }
    else
    {
      ++v_q;
      v_evalvec[v_q] = eval;  seval += eval;
      for (int i = 0; i < ni(); ++i) v_amat.assign(i + 1, v_q, evecv(i));
      if (seval / v_evalvec[0] >= (ldbl)(opt.c())) respca = resok;
      else if (v_q >= qmax) respca = resrmax;
    }
  } while (respca == resnon);
  v_amat.resize(ni() + 1, v_q + 1);
  if (opt.q() < 1 && respca == resrmax) respca = resok;
  return (v_q > 0) ? respca :  reserr;
}

class pcaclass                  // pcaメインクラス
{
private:
  string fname;                 // ファイル名
  idd id;                       // データ識別情報
  ldata datalst;                // データリスト
  void put_resid(ostream& fout);
  // 解析結果のタイトルとデータ識別情報を出力する
  void put_stat(ostream& fout, option& opt, const varclass& var) const;
  // 各種統計量を出力する
  void put_pca(ostream& fout, option& opt, const varclass& var);
  // 主成分計算結果を出力する
  void put_pcadns(option& opt, const varclass& var);
  // 主成分スコアをデータファイル形式で出力する
  void plot_ngpid(ostream& fout);
  // Ngraph用データのタイトルとデータ識別情報を出力する
  void amat2cmat(option& opt, mat<ldbl>& amat, int n, const mat<ldbl>& smat, const vecldbl& evalvec) const;
  // 固有ベクトルを主成分計算用係数に変換する
  void amat2lmat(option& opt, mat<ldbl>& amat, int n, const mat<ldbl>& smat, const vecldbl& evalvec) const;
  // 固有ベクトルを主成分負荷量に変換する
  void plot_plplot(option& opt, const varclass& var);
  // Ngraph用主成分負荷量のプロットデータを出力する
public:
  pcaclass() : fname(), id(), datalst() // コンストラクタ
  {}
  ~pcaclass()                   // デストラクタ
  {}
  int get(const char* c_fname, option& opt); // データを読み込み保存・累積する
  void put(option& opt);                     // 集計結果を出力する
};

void pcaclass::put_resid(ostream& fout)
// 解析結果のタイトルとデータ識別情報を出力する
{
  puttitle(fout,  Title);
  fout << "\n";
  id.put_title(fout);
  id.put_stra(fout);
  for (size_t i = 0; i < id.ni(); ++i)
  { fout << "変数x" << setw(3) << i + 1 << ":" << id.item(i) << "\n"; }
}

void pcaclass::put_stat(ostream& fout, option& opt, const varclass& var) const
// 各種統計量を出力する
{
  putmsdvec(cout, var.n(), var.mevec(), var.vrvec());
  string head = "・計算開始行列:相関行列";
  switch (opt.m())
  {
  case mats : head = "・計算開始行列:積和行列";
    break;
  case matv : head = "・計算開始行列:共分散行列";
    break;
  case matetc : head = "・計算開始行列:単純積和行列";
  }
  if (opt.m() == matetc) putmat(cout, head, mats, var.n(), var.smat());
  else putmat(cout, head, opt.m(), var.n(), var.smat());
}

void pcaclass::put_pca(ostream& fout, option& opt, const varclass& var)
// 主成分計算結果を出力する
{
  constexpr int LENline = 76;   // 区切り線の長さ
  if (opt.res("pc"))
  {
    fout << "\n・主成分計算結果  主成分抽出条件:";
    if (var.resvar() == resrmax) fout << "最大主成分数=" << opt.q();
    else fout << "基準累積寄与率=" << opt.c();
    fout << "\n主成分    固有値  累積固有値      寄与率  累積寄与率      α係数    標準偏差\n";
    putline(fout, LENline);
  }
  mat<ldbl> cmat(var.amat());
  amat2cmat(opt, cmat, var.n(), var.smat(), var.evalvec());
  mat<ldbl> ctmat(cmat);
  for (int i = 0; i < var.ni(); ++i) ctmat.assign(i, 0, 0.0L);
  for (int j = 1; j < var.q(); ++j) ctmat.assign(0, j, 0.0L);
  mat<ldbl> tmat(ctmat.trn() * var.smat());
  mat<ldbl> zsmat(tmat * ctmat.trn());
  int df = var.n() - 1;
  if (opt.res("pc"))
  {
    ldbl slamda = 0.0L;
    for (int k = 1; k <= var.q(); ++k)
    {
      ldbl lamda = var.eval(k);
      ldbl vzz = (opt.m() == matr) ? 1.0L : var.eval(k);
      slamda += lamda;
      if (opt.m() == mats) vzz = lamda / (ldbl)df;
      ldbl sxx = 0.0L;
      for (int i = 1; i <= var.ni(); ++i) sxx += cmat(i, k) * var.sxx(i) * cmat(i, k);
      ldbl alpha = (zsmat(k, k) > 0.0L) ? (ldbl)(var.ni()) / (ldbl)(var.ni() - 1) * (1.0L - sxx / zsmat(k, k)) : 0.0L;
      fout << "z" << setw(3) << k << setw(12) << lamda;
      fout << setw(12) << slamda << setw(12) << lamda / var.eval(0);
      fout << setw(12) << slamda / var.eval(0) << setw(12) << alpha;
      fout << setw(12) << sqrt(vzz) << "\n";
    }
    putline(fout, LENline);
    putcoef(fout, "・固有ベクトル", "主成分", "x", var.amat());
  }
  if (opt.res("pcc")) putcoef(fout, "・主成分スコア計算用係数", "主成分", "x", cmat);
  mat<ldbl> lmat(var.amat());
  amat2lmat(opt, lmat, var.n(), var.smat(), var.evalvec());
  if (opt.res("pcl")) putfl(fout, "・主成分負荷量", "主成分", "x", lmat);
}

void pcaclass::amat2cmat(option& opt, mat<ldbl>& amat, int n, const mat<ldbl>& smat, const vecldbl& evalvec) const
// 固有ベクトルを主成分計算用係数に変換する
{
  int df = n - 1, p = id.ni(), q = amat.c() - 1;
  for (int j = 1; j <= q; ++j)
  {
    ldbl a0 = 0.0, lamda = evalvec[j];
    for (int i = 1; i <= p; ++i)
    {
      ldbl w = 1.0L;            // 共分散行列の時
      switch (opt.m())
      {
      case matr : w = (smat(i, i) > 0.0L) ? (ldbl)df / smat(i, i) : 0.0L;
        break;
      case mats : w = (ldbl)df;
        break;
      case matetc : w = (ldbl)n;
      }
      if (opt.std()) amat.assign(i, j, amat(i, j) * sqrt(w / lamda));
      else amat.assign(i, j, amat(i, j) * sqrt(w));
      a0 -= amat(i, j) * smat(0, i);
    }
    amat.assign(0, j, a0);
  }
}

void pcaclass::amat2lmat(option& opt, mat<ldbl>& amat, int n, const mat<ldbl>& smat, const vecldbl& evalvec) const
// 固有ベクトルを主成分負荷量に変換する
{
  int df = n - 1, p = amat.r() - 1, q = amat.c() - 1;
  for (int i = 1; i <= p; ++i)
  {
    ldbl v = 1.0L;              // 相関行列の場合
    switch (opt.m())
    {
    case matv : v = (df > 0) ? smat(i, i) / (ldbl)df : 0.0L;
      break;
    case matetc :
    case mats : v = smat(i, i);
    }
    if (v <= 0.0L)
    { for (int j = 1; j <= q; ++j) amat.assign(i, j, 0.0); }
    else
    {
      for (int j = 1; j <= q; ++j) amat.assign(i, j, amat(i, j) * sqrt(evalvec[j] / v));
    }
  }
}

void pcaclass::put_pcadns(option& opt, const varclass& var)
// 主成分スコアをデータファイル形式で出力する
{
  ofstream fout;
  if (ofopen(fout, (opt.dnsname()).c_str()) == false) return; // ファイルオープン失敗
  id.put_title(fout, formtag);
  id.put_stra(fout, formtag);
  idd tid(id.title(), var.q());
  ostringstream ostr;
  for (int i = 0; i < var.q(); ++i)
  {
    ostr << "主成分スコア" << i + 1;
    tid.assign_item(i, ostr.str());
    ostr.str("");
  }
  tid.put_item(fout, formtag);
  mat<ldbl> cmat(var.amat());
  amat2cmat(opt, cmat, var.n(), var.smat(), var.evalvec());
  datac score(var.q());
  for (auto& cdatai : datalst)
  {
    if (cdatai.sdata(0) == BLANKs)
    { score = BLANKs; }
    else
    {
      vecdbl fdatavec(id.ni(), 0.0);
      for (int i = 0; i < var.ni(); ++i) fdatavec[i] = cdatai.fdata(i);
      for (int i = 1; i <= var.q(); ++i)
      {
        double z = (double)(cmat(0, i));
        for (int j = 0; j < var.ni(); ++j) z += cmat(j + 1, i) * fdatavec[j];
        score.assign(i - 1, z);
      }
    }
    score.put(fout);
  }
}

void pcaclass::plot_plplot(option& opt, const varclass& var)
// Ngraph用主成分負荷量のプロットデータを出力する
{
  ofstream fout;
  if (ofopen(fout, (opt.ngpname()).c_str()) == false) return; // ファイルオープン失敗
  constexpr int RD = 60;        // 円の分解数
  mat<ldbl> lmat(var.amat());
  amat2lmat(opt, lmat, var.n(), var.smat(), var.evalvec()); // 固有ベクトルを主成分負荷量に変換する
  fout << "-s1 -zx-1,1,0.5 -zy-1,1,0.5\n";
  fout << "%";  puttitle(fout, Title);  fout << "%\n";
  id.put_title(fout, formngp);
  id.put_stra(fout, formngp);
  for (int i = 0; i < var.ni(); ++i)
  { fout << "%変数x" << setw(3) << i + 1 << ":" << id.item(i) << "\n"; }
  fout << "%円を描く時:-dline、プロットを描く時:-dmark,00 -o300\n";
  double dh = 2.0 * M_PI / (double)RD;
  double h = 0.0;
  for (int i = 0; i <= RD; ++i, h += dh) fout << cos(h) << " " <<  sin(h) << "\n";
  for (int j = 1; j < var.q(); ++j)
  {
    for (int k = j + 1; k <= var.q(); ++k)
    {
      fout << NGPdelm << "\n";
      fout << "%x:z" << j << "  y:z" << k << "\n";
      for (int i = 1; i <= var.ni(); ++i) fout << lmat(i, j) << " " << lmat(i, k) << "\n";
    }
  }
}

int pcaclass::get(const char* c_fname, option& opt)
// データを読み込み保存・累積する
{
  if (c_fname == nullptr) fname.erase();
  else fname = c_fname;
  int nc = fgetdata(fname, id, datalst);
  if (nc < 0) throw runtime_error("実行時エラー！");
  if (opt.check(id.ni()) != 0) throw option_error("実行時オプションエラー！");
  return nc;
}

void pcaclass::put(option& opt)
// 集計結果を出力する
{
  varclass var(id.ni());        // 変数オブジェクト
  var.cal(opt, datalst);        // 各種の計算実行
  if (opt.res("res"))
  {
    put_resid(cout);
    if (opt.res("stat")) put_stat(cout, opt, var);
    if (var.resvar() == reserr) cout << "※主成分抽出不可能！\n";
    else put_pca(cout, opt, var);
    cout << "\n" << endl;
  }
  if (var.resvar() == reserr) return;
  if (opt.dnsname() != "") put_pcadns(opt, var);
  if (opt.ngpname() != "") plot_plplot(opt, var);
}

int main(int argc, char* argv[]) // メインルーチン
{
  pcaclass pca;                 // メインクラス変数
  option opt(argv[0]);          // オプション
  try
  {
    opt.setdef();               // オプションの初期設定
    bool fok = true;            // ファイル処理フラッグ
    for (int i = 1; i <= argc; ++i)
    {
      if (i >= argc)
      {
        if (fok) pca.get(nullptr, opt); // 標準入力から読み込む
      }
      else if (*(argv[i]) == CHRopt)
      { opt.set(argv[i] + 1);  fok = true; }
      else
      { pca.get(argv[i], opt);  fok = false; }
    }
    pca.put(opt);
  }
  catch (help_except& helpex)   // ヘルプファイル表示
  { return 0; }
  catch (option_error& opterr)  // オプション指定エラー
  { return 1; }
  catch (runtime_error& rterr)  // 実行時エラー
  {
    disperr(opt.pname(), rterr.what());
    return 2;
  }
  catch (...)                   // 予期せぬ例外
  {
    disperr(opt.pname(), "内容不明なエラー発生！");
    return 255;
  }
  return 0;
}
