// mda.C Version 7.44  Copyleft (c) 杉本典夫
// Time-stamp: <2025-10-31 15:50:38 tac>
//
// 最初の項目を判別項目、以後の項目を変数とした重判別分析(正準分析)を行う。
//
// 使用ライブラリ:libdans.so, libfun.so
// コンパイル法 :g++ -O -s -W mda mda.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 =            // プログラム使用法表示用
{
  "mda  Version 7.44            Copyleft(Vo\\oV) 2025 杉本典夫",
  "書式:mda  [<オプション> …]  [<データファイル> …]",
  "機能:最初の項目を判別項目、以後の項目を変数とした重判別分析をするぜ。",
  "<オプション>",
  "  -dis={euc|mah}:判別用距離の種類指定だぜ。",
  "    euc:ユークリッド平方距離(オプション省略時)  mah:Mahalanobisの汎距離",
  "  -dnsname[=<出力ファイル名>]:DANSデータを出力するファイル名を指定するぜ。",
  "    オプション省略時:DANSデータ出力せず  -dnsname:-dnsname=mda_dns.d",
  "    出力項目:正準スコア",
  "  -f=<ファイル>:指定したファイルからオプションを読み込むぜ。",
  "  -h[=<コマンド名>]:オンラインマニュアルを表示するぜ。(-h:-h=man)",
  "  -res=[-][{def|all|ca|cs|stat|wl}][,…]:解析結果の出力項目を指定するぜ。",
  "    -:-を付けた項目を出力しない(-res=-:全ての項目を出力しない)",
  "    def:オプション省略時と同じ項目を出力する  all:全ての項目を出力する",
  "    ca:正準軸  cs:正準スコア  stat:基礎統計量  wl:WilksのΛ"
};

const string Title = "=== 重判別分析(multiple discriminant analysis, 正準判別分析) ==="; // 表題
const string DNSNAMEdef = "mda_dns.d"; // 省略時DANSファイル名

enum distype                    // 距離の種類
{
  disnon,                       // 指定なし
  diseuc,                       // ユークリッド平方距離
  dismah,                       // Mahalanobisの汎距離
  disetc                        // その他
};

class option                    // オプションクラス
{
private:
  string v_path;                // プログラムのパス
  string v_pname;               // プログラム名
  distype v_dis;                // 距離の種類
  string v_dnsname;             // DANSファイル名
  lstr reslst;                  // 分析結果出力リスト
  mapsb resmap;                 // 分析結果出力の有無マップ
public:
  option(const char* c_pname = nullptr); // コンストラクタ
  ~option()                             // デストラクタ
  {}
  const string& path() const    // プログラムのパスを返す
  { return v_path; }
  const string& pname() const   // プログラム名を返す
  { return v_pname; }
  distype dis() const           // 距離の種類を返す
  { return v_dis; }
  const string& dnsname() const // DANSデータファイル名を返す
  { return v_dnsname; }
  bool res(const string& code); // 結果出力フラッグを返す
  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_dis(diseuc), v_dnsname(), reslst(), resmap()
// コンストラクタ
{
  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] == "dis")          // 距離の種類指定
  {
    if (opt[1].empty() || opt[1] == "euc") v_dis = diseuc;
    else if (opt[1] == "mah") v_dis = dismah;
    else err = badopt(pname(), optstr, USAGE);
  }
  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] == "res")     // 結果の出力指定
  {
    if (strsplit(reslst, opt[1], CHRdelm) < 1) 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 (set_res() < 0) return -3; // -resオプションのエラー
  return 0;                     // エラーなし
}

int option::set_res()           // -resオプションのセット
{
  resmap["res"] = resmap["def"] = resmap["ca"] = resmap["cs"] = resmap["stat"] = resmap["wl"] = true;
  if (reslst.size() > 0)
  {
    resmap["def"] = resmap["ca"] = resmap["cs"] = resmap["stat"] = resmap["wl"] = 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["ca"] = resmap["cs"] = resmap["stat"] = resmap["wl"] = true;
      }
      else if (resitr == resmap.end())
      { return disperr(pname(), "-res=" + resstr + "オプション指定が不適切！"); }
      else
      { resmap[resstr] = true; }
    }
  }
  return 0;                     // エラーなし
}

class gstat                     // 群別統計量クラス
{
private:
  mat<ldbl> v_smat;             // 積和行列
  int v_n;                      // 例数
  vecldbl v_mevec;              // 平均値ベクトル
  vecldbl v_vrvec;              // 積和ベクトル
  mat<ldbl> v_csmat;            // 正準スコア積和行列
  mat<ldbl> v_cimat;            // 正準スコア共分散逆行列
public:
  gstat() : v_smat(), v_n(0), v_mevec(), v_vrvec(), v_csmat(), v_cimat()
  // デフォールトコンストラクタ
  {}
  explicit gstat(size_t p) : v_smat(p + 1), v_n(0), v_mevec(p + 1), v_vrvec(p + 1), v_csmat(1), v_cimat(1)
  // コンストラクタ
  {}
  gstat(const gstat& gstat0) : v_smat(gstat0.v_smat), v_n(gstat0.v_n), v_mevec(gstat0.v_mevec), v_vrvec(gstat0.v_vrvec), v_csmat(gstat0.v_csmat), v_cimat(gstat0.v_cimat)
  // コピーコンストラクタ
  {}
  virtual ~gstat()              // デストラクタ
  {}
  const gstat& operator=(const gstat& gstat0); // コピー代入演算子
  void sigma(const vecldbl& fdatavec);         // データを累積する
  int setsmat();                               // 積和行列をセットする
  void calscore(int p, int g, int q, const mat<ldbl>& amat);
  // 群別正準スコアの積和行列を計算する
  const mat<ldbl>& smat() const // 積和行列を返す
  { return v_smat; }
  ldbl smat(int ridx, int cidx) const // 積和行列の要素を返す
  { return v_smat(ridx, cidx); }
  int n() const                 // 例数を返す
  { return v_n; }
  const vecldbl& mevec() const  // 平均値ベクトルを返す
  { return v_mevec; }
  ldbl me(size_t idx) const     // 平均値を返す
  { return (idx < v_mevec.size()) ? v_mevec[idx] : 0.0L; }
  const vecldbl& vrvec() const // 分散ベクトルを返す
  { return v_vrvec; }
  ldbl vr(size_t idx) const     // 分散を返す
  { return (idx < v_vrvec.size()) ? v_vrvec[idx] : 0.0L; }
  const mat<ldbl>& csmat() const // 正準スコア積和行列を返す
  { return v_csmat; }
  ldbl csmat(int ridx, int cidx) const // 正準スコア積和行列の要素を返す
  { return v_csmat(ridx, cidx); }
  const mat<ldbl>& cimat() const // 正準スコア共分散逆行列を返す
  { return v_cimat; }
  ldbl csm(size_t idx) const    // 正準スコア平均値を返す
  { return (idx < v_csmat.c()) ? v_csmat(0, idx) : 0.0L; }
  ldbl csv(size_t idx) const    // 正準スコア分散を返す
  { return (idx < v_csmat.r() && n() > 1) ? v_csmat(idx, idx) / (ldbl)(n() - 1): 0.0L; }
};

const gstat& gstat::operator=(const gstat& gstat0) // コピー代入演算子
{
  if (&gstat0 == this) return *this; // 自分自身への代入を防ぐ
  v_smat = gstat0.v_smat;  v_n = gstat0.v_n;
  v_mevec = gstat0.v_mevec;  v_vrvec = gstat0.v_vrvec;
  v_csmat = gstat0.v_csmat;  v_csmat = gstat0.v_csmat;
  return *this;
}

void gstat::sigma(const vecldbl& fdatavec) // データを累積する
{
  v_mevec[0] = 1.0L;
  for (size_t i = 0; i < fdatavec.size(); ++i) v_mevec[i + 1] = fdatavec[i];
  for (size_t i = 0; i < v_smat.r(); ++i)
  {
    double x = v_mevec[i];
    for (size_t j = i; j < v_smat.c(); ++j) v_smat.add(i, j, x * v_mevec[j]);
  }
}

int gstat::setsmat()            // 積和行列をセットする
{
  v_n = (int)(v_smat(0, 0));
  for (size_t i = 0; i < v_smat.r() - 1; ++i)
  {
    for (size_t j = i + 1; j < v_smat.c(); ++j) v_smat.assign(j, i, v_smat(i, j));
  }
  v_smat.sweep(0);
  for (size_t i = 0; i < v_smat.r(); ++i) // 平均と分散をセットする
  {
    v_mevec[i] = v_smat(0, i);
    v_vrvec[i] = (v_n > 1) ? v_smat(i, i) / (ldbl)(v_n - 1) : 0.0L;
  }
  return v_n;                   // 例数を返す
}

void gstat::calscore(int p, int g, int q, const mat<ldbl>& amat)
// 群別正準スコアの積和行列を計算する
{
  v_csmat.resize(q + 1, q + 1);
  v_csmat.assign(0, 0, (ldbl)(n())); // 例数
  mat<ldbl> smat(v_smat.sub(1, 1, p, p));
  for (int i = 1; i <= q; ++i)
  {
    ldbl mi = amat(0, i);
    for (int j = 1; j <= p; ++j) mi += v_smat(0, j) * amat(j, i);
    v_csmat.assign(0, i, mi);
    vec<ldbl> aivec(amat.sub(1, i, p, i));
    vec<ldbl> tvec(aivec.trn() * smat);
    for (int j = i; j <= q; ++j)
    {
      vec<ldbl> ajvec(amat.sub(1, j, p, j));
      v_csmat.assign(i, j, inp(tvec, ajvec));
    }
    for (int j = i + 1; j <= q; ++j) v_csmat.assign(j, i, v_csmat(i, j));
  }
  v_cimat = v_csmat.sub(1, 1, q, q);
  if (n() > 1) v_cimat *= 1.0L / (ldbl)(n() - 1);
  v_cimat.inv();                // 逆行列にする
}

using mapsg = map<string, gstat>;         // 群別統計量マップの別名定義
using mapsg_itr = mapsg::iterator;        // 反復子の別名定義
using mapsg_citr = mapsg::const_iterator; // const反復子の別名定義

class varclass                  // 変数クラス
{
private:
  mapsg v_gsmap;                // 群統計量マップ
  mat<ldbl> v_wmat;             // 群内積和行列
  mat<ldbl> v_bmat;             // 群間積和行列
  vecldbl v_mevec;              // 全体の平均値ベクトル
  vecldbl v_vrvec;              // 全体の分散ベクトル
  vecldbl evalvec;              // 固有値ベクトル
  mat<ldbl> v_amat;             // 固有ベクトルの行列
  mat<ldbl> v_samat;            // 累積固有ベクトルの行列
  int v_ni;                     // 項目数
  int v_np;                     // パラメータ数
  int v_nt;                     // 総例数
  int v_qmax;                   // 正順軸最大数
  int v_q;                      // 正順軸数
  restype v_resvar;             // 処理結果の種類
  restype cal_mda();            // 正準軸を抽出する
public:
  varclass(size_t ni) : v_gsmap(), v_wmat(ni), v_bmat(ni), v_mevec(ni, 0.0L), v_vrvec(ni, 0.0L), evalvec(), v_amat(), v_samat(), v_ni(ni), v_np(ni - 1), v_nt(0), v_qmax(0), v_q(0), v_resvar(resnon)
  // コンストラクタ
  {}
  ~varclass()                   // デストラクタ
  {}
  mapsg& gsmap()                // 群統計量マップを返す
  { return v_gsmap; }
  const mat<ldbl>& wmat() const // 群内積和行列を返す
  { return v_wmat; }
  const mat<ldbl>& bmat() const // 群間積和行列を返す
  { return v_bmat; }
  const vecldbl& mevec() const  // 全体の平均値ベクトルを返す
  { return v_mevec; }
  const vecldbl& vrvec() const  // 全体の分散ベクトルを返す
  { return v_vrvec; }
  const mat<ldbl>& amat() const // 固有ベクトルの行列を返す
  { return v_amat; }
  const mat<ldbl>& samat() const // 累積固有ベクトルの行列を返す
  { return v_samat; }
  ldbl tr() const               // トレース(累積固有値)を返す
  { return evalvec[0]; }
  ldbl eval(int idx) const      // 固有値を返す
  { return (0 <= idx && idx <= q()) ? evalvec[idx] : 0.0L; }
  int ni() const                // 項目数を返す
  { return v_ni; }
  int np() const                // パラメータ数を返す
  { return v_np; }
  int g() const                 // 群数を返す
  { return (int)(v_gsmap.size()); }
  int nt() const                // 総例数を返す
  { return v_nt; }
  int qmax() const              // 正順軸最大数を返す
  { return v_qmax; }
  int q() const                 // 正順軸数を返す
  { return v_q; }
  restype resvar() const        // 処理結果の種類を返す
  { return v_resvar; }
  void cal(option& opt, ldata& datalst); // データを集計して各種統計量を求める
  int search_gs(const string& gstr) const; // 群コードをサーチして群インデクスを返す
  int cscore(datac& score, datac& cdata, option& opt) const;
  // 正準スコアを計算し、データとして返す
};

void varclass::cal(option& opt, ldata& datalst)
// データを集計して各種統計量を求める
{
  vecldbl fdatavec(np());
  for (auto& cdatai : datalst)
  {
    bool blank = false;
    for (int j = 0; blank == false && j < np(); ++j)
    {
      fdatavec[j] = cdatai.fdata(j + 1);
      if (fdatavec[j] == BLANKf) blank = true;
    }
    if (cdatai.sdata(0) != BLANKs && blank == false)
    {
      mapsg_itr gitr = v_gsmap.find(cdatai.sdata(0));
      if (gitr == v_gsmap.end()) // 新規群コード
      {
        gstat gstati(np());
        gstati.sigma(fdatavec);
        v_gsmap.insert(make_pair(cdatai.sdata(0), gstati));
      }
      else                      // すでにある群
      { gitr->second.sigma(fdatavec); }
    }
  }
  v_nt = 0;
  for (auto& gstati : v_gsmap)
  {
    for (int j = 0; j < ni(); ++j) v_bmat.add(0, j, gstati.second.smat(0, j));
    v_nt += gstati.second.setsmat();
    v_wmat += gstati.second.smat();
    for (int i = 1; i < ni(); ++i)
    {
      ldbl sigx = gstati.second.n() * gstati.second.me(i);
      for (int j = i; j < ni(); ++j) v_bmat.add(i, j, sigx * gstati.second.me(j));
    }
  }
  v_wmat.assign(0, 0, v_bmat(0, 0));
  for (size_t i = 0; i < v_bmat.r(); ++i)
  {
    for (size_t j = i + 1; j < v_bmat.c(); ++j) v_bmat.assign(j, i, v_bmat(i, j));
  }
  v_bmat.sweep(0);
  v_bmat.assign(0, 0, v_wmat(0, 0));
  for (int k = 1; k < ni(); ++k) v_wmat.assign(0, k, v_bmat(0, k));
  mat<ldbl> stmat(v_wmat + v_bmat);
  for (int i = 0; i <= np(); ++i)
  {
    v_mevec[i] = v_wmat(0, i);
    v_vrvec[i] = stmat(i, i) / (ldbl)(v_nt - 1);
  }
  if (v_gsmap.size() > 1 && v_nt > (int)(v_gsmap.size()) + 1) v_resvar = cal_mda();
  else v_resvar = reserr;
}

restype varclass::cal_mda()     // 正準軸を抽出する
{
  v_qmax = min(g() - 1, np());
  evalvec.resize(v_qmax + 1, 0.0L);
  v_amat.resize(ni(), v_qmax + 1);
  ldbl w = sqrt((ldbl)(nt() - g()));
  mat<ldbl> wmatinv = v_wmat.sub(1, 1, v_wmat.r() - 1, v_wmat.c() - 1);
  wmatinv.inv();
  mat<ldbl> wbmat = wmatinv * v_bmat.sub(1, 1, v_bmat.r() - 1, v_bmat.c() - 1);
  evalvec[0] = 0.0L;
  for (int i = 0; i < np(); ++i) evalvec[0] += wbmat(i, i); // トレースの計算
  v_q = 0;
  vec<ldbl> axv(np()), ayv(np());
  while (v_q < v_qmax && (evalvec[v_q + 1] = eigenud(wbmat, axv, ayv)) > 0.0L)
  {
    ++v_q;
    vec<ldbl> tv = v_wmat.sub(1, 1, v_wmat.r() - 1, v_wmat.c() - 1) * axv;
    ldbl norm = inp(axv, tv);
    axv *= (norm > 0.0L) ? sqrt(1.0L / norm) : 0.0L;
    v_amat.assign(0, v_q, 0.0L);
    for (int i = 1; i <= np(); ++i)
    {
      v_amat.assign(i, v_q, w * axv(i - 1));
      v_amat.add(0, v_q, -v_amat(i, v_q) * v_wmat(0, i));
    }
  }
  v_samat.resize(ni(), v_q + 1);
  for (int i = 1; i <= np(); ++i)
  {
    ldbl sdi = sqrt(v_wmat(i, i) / (ldbl)(nt() - g()));
    for (int j = 1; j <= v_q; ++j) v_samat.assign(i, j, v_amat(i, j) * sdi);
  }
  if (v_q < 1) return reserr;
  return resok;
}

int varclass::search_gs(const string& gstr) const
// 群コードをサーチして群インデクスを返す
{
  mapsg_citr gitr = v_gsmap.begin();
  int gidx = 0;
  bool hit = false;
  while (hit == false && gitr != v_gsmap.end())
  {
    if (gstr == gitr->first) hit = true;
    else { ++gitr;  ++gidx; }
  }
  return hit ? gidx : -1;       // 群インデクスを返す
}

int varclass::cscore(datac& score, datac& cdata, option& opt) const
// 正準スコアを計算し、データとして返す
{
  score.assign(0, cdata.sdata(0));
  vecdbl fdatavec(np(), 0.0);
  bool blank = false;
  for (int j = 0; blank == false && j < np(); ++j)
  {
    fdatavec[j] = cdata.fdata(j + 1);
    if (fdatavec[j] == BLANKf) blank = true;
  }
  if (cdata.sdata(0) == BLANKs || blank) // スコアの計算ができない
  {
    for (size_t i = 1; i < score.size(); ++i) score.assign(i, BLANKs);
    return 0;                   // スコア計算できず
  }
  for (int i = 1; i <= q(); ++i)
  {
    double z = (double)(v_amat(0, i));
    for (int j = 0; j < np(); ++j) z += v_amat(j + 1, i) * fdatavec[j];
    score.assign(i, z);
  }
  int i = q() + 1;
  for (auto& gstati : v_gsmap)
  {
    vec<ldbl> dvec(q());
    for (int j = 1; j <= q(); ++j) dvec.assign(j - 1, score.fdata(j) - gstati.second.csm(j));
    vec<ldbl> dtvec(dvec);
    if (opt.dis() == dismah) dtvec = dvec.trn() * gstati.second.cimat(); // Mahalanobisの汎距離
    score.assign(i, inp(dtvec, dvec));
    ++i;
  }
  return q();                   // 計算できたスコア数を返す
}

class mdaclass                  // mdaメインクラス
{
private:
  string fname;                 // ファイル名
  idd id;                       // データ識別情報
  ldata datalst;                // データリスト
  void put_resid(ostream& fout); // 解析結果のタイトルとデータ識別情報を出力する
  void put_stat(ostream& fout, varclass& var) const;
  // 統計量を出力する
  void put_mdadns(option& opt, varclass& var);
  // 正準スコアをデータファイル形式で出力する
  void put_mdatab(ostream& fout, const varclass& var) const;
  // 正準軸を出力する
  void put_score(ostream& fout, varclass& var) const;
  // 群別正準スコアの基礎統計量を出力する
  void put_dis(ostream& fout, option& opt, const varclass& var);
  // 群別判別率を出力する
  void put_wilks(ostream& fout, const varclass& var) const;
  // WilksのΛ統計量を出力する
public:
  mdaclass() : fname(), id(), datalst() // コンストラクタ
  {}
  ~mdaclass()                   // デストラクタ
  {}
  int get(const char* c_fname, option& opt); // データを読み込み累積する
  void put(option& opt);                     // 集計結果を出力する
};

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

void mdaclass::put_stat(ostream& fout, varclass& var) const
// 統計量を出力する
{
  int gno = 1;
  for (auto& gstati : var.gsmap())
  {
    fout << "\n○群" << gno << ":" << id.item(0) << "=" << gstati.first;
    putmsdvec(fout, gstati.second.n(), gstati.second.mevec(), gstati.second.vrvec());
    ++gno;
  }
  fout << "\n○全体";
  putmsdvec(fout, var.nt(), var.mevec(), var.vrvec());
  putmat(fout, "・群内相関行列(Within)", matr, var.nt(), var.wmat());
}

void mdaclass::put_mdadns(option& opt, 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(), 1 + var.q() + var.g());
  tid.assign_item(0, id.item(0));
  ostringstream ostr;
  for (int i = 1; i <= var.q(); ++i)
  {
    ostr << "正準スコア" << i;
    tid.assign_item(i, ostr.str());
    ostr.str("");
  }
  int i = var.q() + 1, gno = 1;
  for (auto& gstati : var.gsmap())
  {
    ostr << "群" << gno;
    if (opt.dis() == diseuc) ostr << "の重心からのユークリッド平方距離";
    else ostr << "の重心からのMahalanobisの汎距離";
    tid.assign_item(i, ostr.str());
    ostr.str("");
    gstati.second.calscore(var.np(), var.g(), var.q(), var.amat());
    ++i;  ++gno;
  }
  tid.put_item(fout, formtag);
  datac score(1 + var.q() + var.g());
  for (auto& cdatai : datalst)
  {
    var.cscore(score, cdatai, opt);
    score.put(fout);
  }
}

void mdaclass::put_mdatab(ostream& fout, const varclass& var) const
// 正準軸を出力する
{
  fout << "\n・重判別分析計算結果\n";
  fout << "正準軸          固有値        累積固有値            寄与率        累積寄与率\n";
  constexpr int LENline = 76;   // 区切り線の長さ
  putline(fout, LENline);
  ldbl slam = 0.0L, tr = var.tr();
  for (int i = 1; i <= var.q(); ++i)
  {
    ldbl lam = var.eval(i);
    slam += lam;
    fout << "z" << setw(3) << i << setw(18) << lam << setw(18) << slam;
    fout << setw(18) << lam / tr << setw(18) << slam / tr << "\n";
  }
  putline(fout, LENline);
}

void mdaclass::put_score(ostream& fout, varclass& var) const
// 群別正準スコアの基礎統計量を出力する
{
  constexpr int LENline = 72;   // 区切り線の長さ
  int df = var.nt() - var.g();
  for (int i = 1; i <= var.q(); ++i)
  {
    fout << "\n・正準スコアz" << i << "\n";
    putline(fout, LENline);
    ldbl ss = (ldbl)df;
    int gno = 1;
    for (auto& gstati : var.gsmap())
    {
      gstati.second.calscore(var.np(), var.g(), var.q(), var.amat());
      ldbl m = gstati.second.csm(i), v = gstati.second.csv(i);
      fout << "群" << setw(2) << gno << ":";
      putmsd(fout, gstati.second.n(), m, v);
      ss += gstati.second.n() * m * m;  ++gno;
    }
    putline(fout, LENline);
    fout << "全体:";
    putmsd(fout, var.nt(), 0.0, ss / (ldbl)(var.nt() - 1));
  }
}

void mdaclass::put_dis(ostream& fout, option& opt, const varclass& var)
// 判別結果を出力する
{
  datac score(1 + var.q() + var.g());
  ftable ftab(var.g(), var.g());
  for (auto& cdatai : datalst)
  {
    int gidx = var.search_gs(cdatai.sdata(0));
    if (gidx >= 0 && var.cscore(score, cdatai, opt) > 0)
    {
      double dmin = FDmax;
      int midx = 0;
      for (int i = 1 + var.q(); i < 1 + var.q() + var.g(); ++i)
      {
        if (score.fdata(i) < dmin)
        { dmin = score.fdata(i);  midx = i - 1 - var.q(); }
      }
      ftab.f(gidx + 1, midx + 1, ftab.f(gidx + 1, midx + 1) + 1);
    }
  }
  ostringstream ostr;
  string tstr("群＼判別");
  int npad = tstr.size() - strclen(tstr, CHARwide);
  int wr = 8, wc = 5, tn = 0;
  ostr << left << setw(wr + npad) << tstr << right;
  ftab.hed(ostr.str());
  for (int i = 1; i <= var.g(); ++i)
  {
    ostr.str("");
    tstr = "群";
    npad = tstr.size() - strclen(tstr, CHARwide);
    ostr << setw(wc + npad) << tstr << i;
    ftab.ccode(i, ostr.str());
    ostr.str("");
    ostr << tstr << i;
    ftab.rcode(i, ostr.str());
    tn += ftab.f(i, i);
  }
  ftab.set();
  fout << "\n・正準スコアによる判別結果";
  if (opt.dis() == diseuc) fout << "(ユークリッド平方距離が最短の群)";
  else fout << "(Mahalanobisの汎距離が最短の群)";
  putftable(fout, ftab, tabfprr);
  fout << "全体の判別率=" << (double)tn / (double)(ftab.ftotal()) * 100.0 << "%(" << tn << "/" << ftab.ftotal() << ")\n";
}

void mdaclass::put_wilks(ostream& fout, const varclass& var) const
// WilksのΛ統計量を出力する
{
  ldbl lamd = 1.0L, f = 0.0L;
  for (int i = 1; i <= var.q(); ++i) lamd *= (var.eval(i) + 1.0L);
  lamd = 1.0 / lamd;
  int df1 = 0, df2 = 0, w = var.nt() - var.g();
  if (var.g() == 2)
  {
    df1 = var.np();  df2 = w - var.np() + 1;
    f = (1.0 - lamd) / lamd * (w - var.np() + 1) / var.np();
  }
  else if (var.g() == 3)
  {
    df1 = 2 * var.np();  df2 = 2 * (w - var.np() + 1);
    f = (1.0 - sqrt(lamd)) / sqrt(lamd) * (w - var.np() + 1) / var.np();
  }
  else if (var.np() == 2 && var.g() > 3)
  {
    df1 = 2 * (var.g() - 1);  df2 = 2 * (w - 1);
    f = (1.0L - sqrt(lamd)) / sqrt(lamd) * (ldbl)(w - 1) / (ldbl)(var.g() - 1);
  }
  else
  {
    df1 = (var.g() - 1) * var.np();  df2 = 0;
    f = -((ldbl)w - (ldbl)(var.np() - var.g()) * 0.5L) * log(lamd);
  }
  fout << "\n・WilksのΛ(ラムダ統計量)=" << lamd << "  相関比:η=" << sqrt(1.0L - lamd) << "\n";
  if (df2 > 0)
  {
    ldbl pr = ftop(f, df1, df2);
    fout << "  F=" << f << "  第1自由度=" << df1 << "  第2自由度=" << df2 << "  有意確率:p=" << pr << pind(pr) << endl;
  }
  else
  {
    ldbl pr = ctop(f, df1);
    fout << "  χ2乗=" << f << "  自由度=" << df1 << "  有意確率p=" << pr << pind(pr) << endl;
  }
}

int mdaclass::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 mdaclass::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, var);
    if (var.resvar() == reserr)
    { cout << "\n※重判別関数計算不可能！\n\n" << endl; }
  }
  if (var.resvar() == reserr) return;
  if (opt.res("res"))
  {
    if (opt.res("ca"))
    {
      put_mdatab(cout, var);
      putcoef(cout, "・正準係数(正準スコア計算用)", "正準軸", "x", var.amat());
      putcoef(cout, "・標準正準係数(群内分散で標準化)", "正準軸", "x", var.samat());
    }
    if (opt.res("cs"))
    {
      put_score(cout, var);
      put_dis(cout, opt, var);
    }
    if (opt.res("wl")) put_wilks(cout, var);
    cout << "\n" << endl;
  }
  if (opt.dnsname() != "") put_mdadns(opt, var);
}

int main(int argc, char* argv[]) // メインルーチン
{
  mdaclass mda;                 // メインクラス変数
  option opt(argv[0]);          // オプション
  try
  {
    opt.setdef();               // オプションの初期設定
    bool fok = true;            // ファイル処理フラッグ
    for (int i = 1; i <= argc; ++i)
    {
      if (i >= argc)
      {
        if (fok) mda.get(nullptr, opt); // 標準入力から読み込む
      }
      else if (*(argv[i]) == CHRopt)
      { opt.set(argv[i] + 1);  fok = true; }
      else
      { mda.get(argv[i], opt);  fok = false; }
    }
    mda.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;
}
