// ancova2.C Version 7.44  Copyleft (c) 杉本典夫
// Time-stamp: <2025-10-30 11:10:24 tac>
//
// 最初の項目を群項目、以後の項目を共変数Xと集計項目Yにし、
// 時期をダミー変数化した共変数とする二元配置型共分散分析
//
// 使用ライブラリ:libdans.so, libfun.so
// コンパイル法  :g++ -O -s -W ancova2.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 =            // プログラム使用法表示用
{
  "ancova2  Version 7.44        Copyleft(Vo\\oV) 2025 杉本典夫",
  "書式:ancova2 [<オプション> …]  [<データファイル> …]",
  "機能:最初の項目を群項目、以後の項目を共変数Xと集計項目Yにし、",
  "     時期を共変数にした二元配置型共分散分析をするぜ。",
  "<オプション>",
  "  -cc[=<信頼係数>]:信頼係数を指定するぜ。(0≦信頼係数<1)",
  "    -cc:-cc=0.95(オプション省略時)",
  "  -f=<ファイル>:指定したファイルからオプションを読み込むぜ。",
  "  -h[=<コマンド名>]:オンラインマニュアルを表示するぜ。(-h:-h=man)",
  "  -nx=<共変数の個数>:共変数の個数を指定するぜ。",
  "  -res=[-][{def|all|atab|bt|gstat|m[{a|d}]{b|f|s|t}|pd|tstat}]",
  "       [,…]:解析結果の出力項目を指定する。",
  "    -:-を付けた項目を出力しない(-res=-:全ての項目を出力しない)",
  "    def:オプション省略時と同じ項目を出力する  all:全ての項目を出力する",
  "    atab:共分散分析表  bt:回帰係数の検定  gstat:群別統計量",
  "    m[{a|d}]{b|f|s|t}:修正群差と非平行性の多重比較",
  "    pd:非平行性を残差にプールした分散分析表  tstat:全体の統計量",
  "    オプション省略時:自動判定",
  "    多重比較 b:Bonferroni型  f:Fisher型  s:Scheffe型  t:Tukey型",
  "  -t=<時間1>,<時間2>[,<時間3>,…]:時期項目の値を指定するぜ。",
  "    オプション省略時:時期はダミー変数"
};

const string Title = "=== 二元配置型共分散分析(ANCOVA2) ==="; // 表題

class option                    // オプションクラス
{
private:
  string v_path;                // プログラムのパス
  string v_pname;               // プログラム名
  double v_cc;                  // 信頼係数
  mctype v_ma;                  // 修正群差の多重比較法
  mctype v_md;                  // 非平行性の多重比較法
  lstr reslst;                  // 解析結果出力リスト
  mapsb resmap;                 // 解析結果出力の有無マップ
  int v_nx;                     // 共変数の個数
  lstr tlst;                    // 時期変数の値リスト
  vecdbl v_tvec;                // 時期変数の値ベクトル
public:
  option(const char* c_pname = nullptr); // コンストラクタ
  ~option()                              // デストラクタ
  {}
  const string& path() const    // プログラムのパスを返す
  { return v_path; }
  const string& pname() const   // プログラム名を返す
  { return v_pname; }
  double cc() const             // 信頼係数を返す
  { return v_cc; }
  mctype ma() const             // 修正群差の多重比較法を返す
  { return v_ma; }
  mctype md() const             // 非平行性の多重比較法を返す
  { return v_md; }
  bool res(const string& code); // 結果出力フラッグを返す
  int nx() const                // 共変数の個数を返す
  { return v_nx; }
  const vecdbl& tvec() const    // 時期変数の値ベクトルを返す
  { return v_tvec; }
  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_cc(CCdef), v_ma(mctuk), v_md(mcbon), reslst(), resmap(), v_nx(0), tlst(), v_tvec()
// コンストラクタ
{
  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] == "cc")           // 信頼係数の指定
  {
    v_cc = opt[1].empty() ? CCdef : strtodf(opt[1]);
    if (v_cc < 0.0 || 1.0 <= v_cc) err = badopt(pname(), optstr, USAGE, 0.0, 1.0);
  }
  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] == "nx")      // 共変数の個数指定
  {
    char* ptr = nullptr;
    v_nx = strtol(opt[1].c_str(), &ptr, 0);
    if (v_nx < 0) err = badopt(pname(), optstr, USAGE, 0.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] == "t")       // 時期変数の値指定
  {
    if (strsplit(tlst, 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)    // オプションのチェックとセット
{
  int p = ni - 1 - nx();
  if (p < 1) throw runtime_error("項目が足らんぜ！");
  if (set_res() < 0) return -3; // -resオプションのエラー
  if (tlst.size() < 1) return 0; // 時期はダミー変数
  v_tvec.resize(p);
  if (settvec(v_tvec, tlst, p) < p) return disperr(pname(), "項目時間が間違ってるぜ！");
  return 0;                     // エラーなし
}

int option::set_res()           // -resオプションのセット
{
  resmap["res"] = resmap["def"] = resmap["atab"] = resmap["gstat"] = resmap["tstat"] = true;
  resmap["bt"] = resmap["ma"] = resmap["md"] = resmap["pd"] = false;
  if (reslst.size() > 0)
  {
    resmap["def"] = resmap["atab"] = resmap["gstat"] = resmap["tstat"] = 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") // 省略時出力
      {
        resmap["def"] = resmap["atab"] = resmap["gstat"] = resmap["tstat"] = true;
      }
      else if (resstr == "all") // 全ての結果を出力する
      {
        resmap["def"] = resmap["atab"] = resmap["gstat"] = resmap["tstat"] = resmap["bt"] = resmap["ma"] = resmap["md"] = resmap["pd"] = true;
      }
      else if (resstr == "m")      // -res=mと指定された時
      { resmap["ma"] = resmap["md"] = true; }
      else if (resstr[0] == 'm' && resstr[1] == 'a' && resstr.size() > 2)
      {
        switch (resstr[2])
        {
        case 'b' : v_ma = mcbon;
          break;
        case 'f' : v_ma = mcfis;
          break;
        case 's' : v_ma = mcsch;
          break;
        case 't' : v_ma = mctuk;
          break;
        default : return disperr(pname(), "-res=" + resstr + "オプション指定が不適切！");
        }
        resmap["ma"] = true;
      }
      else if (resstr[0] == 'm' && resstr[1] == 'd' && resstr.size() > 2)
      {
        switch (resstr[2])
        {
        case 'b' : v_md = mcbon;
          break;
        case 'f' : v_md = mcfis;
          break;
        case 's' : v_md = mcsch;
          break;
        default : return disperr(pname(), "-res=" + resstr + "オプション指定が不適切！");
        }
        resmap["md"] = true;
      }
      else if (resitr == resmap.end())
      { return disperr(pname(), "-res=" + resstr + "オプション指定が不適切！"); }
      else
      { resmap[resstr] = true; }
    }
  }
  return 0;                     // エラーなし
}

class gstat                     // 群別統計量クラス
{
private:
  option* optptr;               // オプションへのポインター
  int v_p;                      // 時期数
  int v_dim;                    // ベクトルと行列の次数(x+時期変数の個数+2)
  vector<statv<double>> xstatvec; // xの基礎統計量ベクトル
  vector<statv<double>> ystatvec; // yの基礎統計量ベクトル
  mat<double> smat;             // 積和行列
  vecdbl mevec;                 // xとyの平均ベクトル
  vecdbl ssvec;                 // xとyの平方和ベクトル
public:
  gstat() : optptr(0), v_p(0), v_dim(0), xstatvec(), ystatvec(), smat(), mevec(), ssvec()
  // デフォールトコンストラクタ
  {}
  gstat(option& opt0, int p0); // コンストラクタ
  gstat(const gstat& gstat0) : optptr(gstat0.optptr), v_p(gstat0.v_p), v_dim(gstat0.v_dim), xstatvec(gstat0.xstatvec), ystatvec(gstat0.ystatvec), smat(gstat0.smat), mevec(gstat0.mevec), ssvec(gstat0.ssvec)
  // コピーコンストラクタ
  {}
  virtual ~gstat()              // デストラクタ
  {}
  const gstat& operator=(const gstat& gstat0); // コピー代入演算子
  const gstat& operator+=(const gstat& gstat0); // 加算代入演算子
  int sigma(const vecdbl& fdatavec); // データを累積してデータ数を返す
  bool setsmat(bool sets); // 平均と積和行列をセットして計算可能性を返す
  void calreg();           // 回帰式を計算する
  int p() const            // 時期数を返す
  { return v_p; }
  int nx() const                // 共変数xの個数を返す
  { return optptr->nx(); }
  int dim() const    // ベクトルと行列の次数を返す(x+時期変数の個数+2)
  { return v_dim; }
  const vecdbl& tvec() const    // 時期変数の値ベクトルを返す
  { return optptr->tvec(); }
  const statv<double>& xstat(int idx) const // xの基礎統計量を返す
  { return xstatvec[idx]; }
  const statv<double>& ystat(int idx) const // yの基礎統計量を返す
  { return ystatvec[idx]; }
  int n() const                 // 全例数を返す
  { return ystatvec[p()].n(); }
  double me(int idx) const      // 平均値を返す
  { return (0 <= idx && idx < dim()) ? mevec[idx] : 0.0; }
  double my() const             // yの平均値を返す
  { return mevec[dim() - 1]; }
  double ss(int idx) const      // 平方和を返す
  { return (0 <= idx && idx < dim()) ? ssvec[idx] : 0.0; }
  double syy() const            // yの平方和を返す
  { return ssvec[dim() - 1]; }
  double vy() const             // yの全分散を返す
  { return ystatvec[p()].vr(); }
  double b(int idx) const       // 偏回帰係数を返す
  { return (0 <= idx && idx < dim() - 1) ? smat(idx, dim() - 1) : 0.0; }
  double vbinv(int idx) const   // 回帰係数の分散要素を返す
  { return (0 <= idx && idx < dim() - 1) ? smat(idx, idx) : 0.0; }
  double sr() const             // 回帰残差を返す
  { return smat(dim() - 1, dim() - 1); }
  double vxr(const vec<double>& dxv) const; // 多重比較用dV^(-1)dを返す
};

gstat::gstat(option& opt0, int p0) : optptr(&opt0), v_p(p0)
// コンストラクタ
{
  if (optptr == 0) return;           // オプションポインターなし
  int dxp = optptr->nx() + p0 + 1; // 時期がダミー変数の時の次数
  v_dim = ((optptr->tvec()).size() > 0) ? optptr->nx() + 3 : dxp;
  xstatvec.resize(optptr->nx());  ystatvec.resize(p0 + 1);
  smat.resize(dxp, dxp);  mevec.resize(dxp, 0.0);  ssvec.resize(dxp, 0.0);
}

const gstat& gstat::operator=(const gstat& gstat0) // コピー代入演算子
{
  if (&gstat0 == this) return *this; // 自分自身への代入を防ぐ
  optptr = gstat0.optptr;  v_p = gstat0.v_p;  v_dim = gstat0.v_dim;
  xstatvec = gstat0.xstatvec;  ystatvec = gstat0.ystatvec;
  smat = gstat0.smat;  mevec = gstat0.mevec;  ssvec = gstat0.ssvec;
  return *this;                 // 自分自身を返す
}

const gstat& gstat::operator+=(const gstat& gstat0) // 加算代入演算子
{
  for (int j = 0; j < optptr->nx(); ++j) xstatvec[j] += gstat0.xstatvec[j];
  for (int j = 0; j < p(); ++j) ystatvec[j] += gstat0.ystatvec[j];
  smat += gstat0.smat;
  return *this;                 // 自分自身を返す
}

int gstat::sigma(const vecdbl& fdatavec) // データを累積してデータ数を返す
{
  mevec[0] = 1.0;
  for (int i = 0; i < optptr->nx(); ++i)
  {
    xstatvec[i].sigma(fdatavec[i]);
    mevec[1 + i] = fdatavec[i];
  }
  int nd = 0;
  for (int j = 0; j < p(); ++j)
  {
    int i = optptr->nx() + j;
    if (fdatavec[i] != BLANKf)
    {
      ystatvec[j].sigma(fdatavec[i]);  ++nd;
      if ((optptr->tvec()).size() > 0) // 時期は連続変数
      { mevec[1 + optptr->nx()] = (optptr->tvec())[j]; }
      else
      {
        for (int k = 1; k < p(); ++k)
        {
          if (optptr->nx() + k == i) mevec[optptr->nx() + k] = 1.0;
          else mevec[optptr->nx() + k] = 0.0;
        }
      }
      mevec[dim() - 1] = fdatavec[i];
      for (int k = 0; k < dim(); ++k)
      {
        double xk = mevec[k];
        for (int l = k; l < dim(); ++l) smat.add(k, l, xk * mevec[l]);
      }
    }
  }
  return nd;                    // 累積したyのデータ数を返す
}

bool gstat::setsmat(bool sets) // 平均と積和行列をセットして例数を返す
{
  for (int i = 0; i < p(); ++i) ystatvec[p()] += ystatvec[i];
  if (sets)                     // 単純積和行列を積和行列にする
  {
    for (int i = 0; i < dim(); ++i)
    { for (int j = i + 1; j < dim(); ++j) smat.assign(j, i, smat(i, j)); }
    smat.sweep(0);
  }
  for (int i = 0; i < dim(); ++i)
  { mevec[i] = smat(0, i);  ssvec[i] = smat(i, i); }
  return (ystatvec[p()].n() < dim() - 1) ? false : true; // 計算可能性を返す
}

void gstat::calreg()            // 回帰式を計算する
{
  for (int i = 1; i < dim() - 1; ++i) smat.sweep(i);
}

double gstat::vxr(const vec<double>& dxv) const // 多重比較用dV^(-1)dを返す
{
  int dd = dxv.dim();
  vec tv = smat.sub(1, 1, dd, dd) * dxv;
  return inp(dxv, tv);
}

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

class varclass                 // 変数クラス
{
private:
  int v_ni;                    // 項目数
  mapsg v_gsmap;               // 群統計量マップ
  gstat v_gstatall;            // 全体統計量オブジェクト
  gstat v_gstatcom;            // 共通統計量オブジェクト
  double v_sbc;                // 共通回帰平方和
  double v_sa;                 // 群差平方和
  double v_sbt;                // 全体回帰平方和
  double v_saa;                // 修正群差平方和
  double v_sd;                 // 非平行性平方和
  double v_sr;                 // 残差平方和
  void cal_ancova2();          // 共分散分析用各種パラメターを計算する
public:
  varclass(option& opt0, size_t ni); // コンストラクタ
  ~varclass()                   // デストラクタ
  {}
  int ni() const                // 項目数を返す
  { return v_ni; }
  int nt() const                // 全例数を返す
  { return v_gstatall.n(); }
  size_t ng() const             // 群数を返す
  { return v_gsmap.size(); }
  const mapsg& gsmap() const    // 群統計量マップを返す
  { return v_gsmap; }
  const gstat& gstatall() const // 全体統計量オブジェクトを返す
  { return v_gstatall; }
  const gstat& gstatcom() const // 共通統計量オブジェクトを返す
  { return v_gstatcom; }
  int dfb() const               // 回帰自由度を返す
  { return v_gstatall.dim() - 2; }
  int dfa() const               // 群自由度を返す
  { return ng() - 1; }
  int dfd() const               // 非平行性自由度を返す
  { return dfb() * dfa(); }
  int dfr() const               // 残差自由度を返す
  { return dft() - dfa() - ng() * dfb(); }
  int dft() const               // 全体自由度を返す
  { return v_gstatall.n() - 1; }
  double sbc() const            // 共通回帰平方和を返す
  { return v_sbc; }
  double sa() const             // 群差平方和を返す
  { return v_sa; }
  double sbt() const            // 全体回帰平方和を返す
  { return v_sbt; }
  double saa() const            // 修正群差平方和を返す
  { return v_saa; }
  double sd() const             // 非平行性平方和を返す
  { return v_sd; }
  double sr() const             // 残差平方和を返す
  { return v_sr; }
  double st() const             // 全体平方和を返す
  { return v_gstatall.syy(); }
  void cal(option& opt, ldata& datalst);
  // データを集計して各種統計量を求める
};

varclass::varclass(option& opt0, size_t ni0) : v_ni(ni0), v_gsmap(), v_gstatall(), v_gstatcom(), v_sbc(0.0), v_sa(0.0), v_sbt(0.0), v_saa(0.0), v_sd(0.0), v_sr(0.0)
// コンストラクタ
{
  int p = v_ni - 1 - opt0.nx();
  gstat tgstat(opt0, p);        // 一時的統計量オブジェクト
  v_gstatall = tgstat;  v_gstatcom = tgstat;
}

void varclass::cal_ancova2()   // 共分散分析用各種パラメターを計算する
{
  v_sa = 0.0;
  double ssb = 0.0;
  for (auto& gstati : v_gsmap)
  {
    v_sa += gstati.second.n() * gstati.second.my() * gstati.second.my();
    ssb += gstati.second.syy() - gstati.second.sr();
  }
  v_sa -= v_gstatall.n() * v_gstatall.my() * v_gstatall.my();
  if (v_sa < FDeps) v_sa = 0.0;
  v_sbt = st() - v_gstatall.sr();
  if (v_sbt < FDeps) v_sbt = 0.0;
  v_sbc = v_gstatcom.syy() - v_gstatcom.sr();
  if (v_sbc < FDeps) v_sbc = 0.0;
  v_saa = v_sa + v_sbc - v_sbt;
  if (v_saa < FDeps) v_saa = 0.0;
  v_sd = ssb - v_sbc;
  if (v_sd < FDeps) v_sd = 0.0;
  v_sr = st() - v_sa - ssb;
  if (dfr() < 1 || v_sr < FDeps) v_sr = 0.0;
}

void varclass::cal(option& opt, ldata& datalst)
// データを集計して各種統計量を求める
{
  int nxp = ni() - 1, p = ni() - 1 - opt.nx();
  vecdbl fdatavec(nxp);
  for (auto& cdatai : datalst)
  {
    bool xok = true;
    for (int j = 0; xok && j < opt.nx(); ++j) // 共変数のチェック
    {
      if (cdatai.fdata(j + 1) == BLANKf) xok = false;
      else fdatavec[j] = cdatai.fdata(j + 1);
    }
    if (cdatai.sdata(0) != BLANKs && xok)
    {
      for (int j = opt.nx(); j < nxp; ++j) fdatavec[j] = cdatai.fdata(j + 1);
      mapsg_itr gitr = v_gsmap.find(cdatai.sdata(0));
      if (gitr == v_gsmap.end()) // 新規群コード
      {
        gstat gstati(opt, p);
        gstati.sigma(fdatavec);
        v_gsmap.insert(make_pair(cdatai.sdata(0), gstati));
      }
      else                      // すでにある群
      { gitr->second.sigma(fdatavec); }
    }
  }
  for (auto& gstati : v_gsmap)
  {
    v_gstatall += gstati.second;
    bool setsmat = gstati.second.setsmat(true);
    v_gstatcom += gstati.second;
    if (setsmat) gstati.second.calreg();
  }
  v_gstatall.setsmat(true);  v_gstatcom.setsmat(false);
  v_gstatall.calreg();  v_gstatcom.calreg();
  cal_ancova2();
}

class ancova2class              // ancova2メインクラス
{
private:
  option opt;                   // オプション
  string fname;                 // ファイル名
  idd id;                       // データ識別情報
  ldata datalst;                // データリスト
  void put_resid(ostream& fout);
  // 解析結果のタイトルとデータ識別情報を出力する
  void put_mreg(ostream& fout, const varclass& var);
  // 群別基礎統計量と回帰式を出力する
  void put_mregtab(ostream& fout, const gstat& gstati, const gstat& gstatcom);
  // 基礎統計量と回帰式の表を出力する
  void put_bt(ostream& fout, const varclass& var);
  // void put_bt(ostream& fout, gstat& gstatcom, int dfr, double sr);
  // 回帰係数の検定結果を出力する
  void put_ancova2tab(ostream& fout,double sbc, int dfb, double sa, int dfa, double sbt, double saa, double sd, int dfd, double sr, int dfr, double st, int dft) const;
  // 共分散分析表を出力する
  void put_macmp(ostream& fout, double sr, int dfr, const varclass& var);
  // 修正群差の多重比較を出力する
  void put_maci(ostream& fout, double sr, int dfr, const varclass& var);
  // 修正群差の信頼区間を出力する
  void put_mdcmp(ostream& fout, double sr, int dfr, const varclass& var);
  // 非平行性の多重比較を出力する
public:
  ancova2class(const char* c_pname = nullptr) : opt(c_pname), fname(), id(), datalst()
  // コンストラクタ
  {}
  ~ancova2class()               // デストラクタ
  {}
  const string& pname() const   // プログラム名を返す
  { return opt.pname(); }
  void init()      // 初期設定を行う(コンストラクタ中ではできない処理)
  { opt.setdef(); }
  void setopt(const string& optstr, char devc = OPTdev) // オプションを設定する
  { opt.set(optstr, devc); }
  int get(const char* c_fname); // データを読み込み累積する
  void put(int nc);             // 集計結果を出力する
};

void ancova2class::put_resid(ostream& fout)
// 解析結果のタイトルとデータ識別情報を出力する
{
  puttitle(fout, Title);
  fout << "\n";
  id.put_title(fout);
  id.put_stra(fout);
  fout << "群項目    :" << id.item(0) << "\n";
  for (int i = 1; i < int(id.ni()); ++i)
  {
    if (i <= opt.nx()) fout << "共変数" << setw(4) << i;
    else fout << "集計項目" << setw(2) << i - opt.nx();
    fout << ":" << id.item(i) << "\n";
  }
}

void ancova2class::put_mreg(ostream& fout, const varclass& var)
// 群別基礎統計量と回帰式を出力する
{
  int nt = (var.gstatall()).n(), ng = (var.gsmap()).size();
  int p = id.ni() - 1 - opt.nx();
  if (opt.res("gstat"))
  {
    int gno = 1;
    for (auto& gstati : var.gsmap())
    {
      fout << "\n・群" << setw(5) << gno;
      fout << ":" << id.item(0) << "=" << gstati.first << "\n";
      put_mregtab(fout, gstati.second, var.gstatcom());
      ++gno;
    }
  }
  if (opt.res("tstat"))
  {
    fout << "\n・全体\n";
    put_mregtab(fout, var.gstatall(), var.gstatcom());
  }
}

void ancova2class::put_mregtab(ostream& fout, const gstat& gstati, const gstat& gstatcom)
// 基礎統計量と回帰式の表を出力する
{
  constexpr int LENline = 77;   // 区切り線の長さ
  constexpr int Nterm = 4;      // 1行に出力する項数
  putline(fout, LENline);
  for (int i = 0; i < gstati.nx(); ++i)
  {
    fout << "x" << left << setw(8) << i + 1 << right << ":";
    putmsd(fout, gstati.xstat(i).n(), gstati.xstat(i).mean(), gstati.xstat(i).vr());
  }
  for (int j = 0; j < gstati.p(); ++j)
  {
    fout << "y:t" << left << setw(6);
    if ((gstati.tvec()).size() > 0) fout << (gstati.tvec())[j];
    else fout << j;
    fout << right << ":";
    putmsd(fout, gstati.ystat(j).n(), gstati.ystat(j).mean(), gstati.ystat(j).vr());
  }
  putline(fout, LENline);
  fout << "y        :";
  putmsd(fout, gstati.n(), gstati.my(), gstati.vy());
  double bc0 = gstati.my(), syy = gstati.syy();
  double sb = syy - gstati.sr();
  double r2 = (syy > 0.0) ? sb / syy : 0.0;
  fout << "群別回帰式:y=" << setw(13) <<  gstati.b(0);
  for (int j = 1; j < gstati.dim() - 1; ++j)
  {
    if ((j + 1) % Nterm == 1) fout << "\n" << setw(13) << "";
    fout << showpos << setw(9) << gstati.b(j) << noshowpos;
    if (j <= gstati.nx())
    { fout << "x" << left << setw(2) << j << right; }
    else
    {
      if ((gstati.tvec()).size() > 0) fout << "t";
      else fout << "t" << left << setw(2) << j - gstati.nx();
      fout << right;
    }
    bc0 -= gstatcom.b(j) * gstati.me(j);
  }
  fout << "\n共通回帰式:y=" << setw(13) << bc0;
  for (int j = 1; j < gstatcom.dim() - 1; ++j)
  {
    if ((j + 1) % Nterm == 1) fout << "\n" << setw(13) << "";
    fout << showpos << setw(9) << gstatcom.b(j) << noshowpos;
    if (j <= gstatcom.nx())
    { fout << "x" << left << setw(2) << j << right; }
    else
    {
      if ((gstatcom.tvec()).size() > 0) fout << "t";
      else fout << "t" << left << setw(2) << j - gstatcom.nx();
      fout << right;
    }
  }
  double r = sqrt(r2);
  int dfb = gstati.dim() - 2;
  int dfr = gstati.n() - 1 - dfb;
  double f = (gstati.sr() > 0.0) ? sb / (double)dfb / gstati.sr() * (double)dfr : 0.0;
  ldbl pr = ftop(f, dfb, dfr);
  fout << "\n群別回帰式の寄与率:r2乗=" << r2 << "  r=" << r << "  有意確率:p=" << pr << pind(pr) << endl;
}

void ancova2class::put_bt(ostream& fout, const varclass& var)
// 回帰係数の検定結果を出力する
{
  constexpr int LENline = 78;   // 区切り線の長さ
  int p = id.ni() - 2;
  double vr = (var.dfr() > 0) ? var.sr() / var.dfr() : 0.0;
  fout << "\n";
  fout << "・共通回帰式の偏回帰係数\n";
  fout << "変数  偏回帰係数   標準誤差" << setw(5) << opt.cc() * 100.0 << "%CI下限        上限         F値 有意確率p値\n";
  putline(fout, LENline);
  for (int j = 1; j < (var.gstatcom()).dim() - 1; ++j)
  {
    double bj = (var.gstatcom()).b(j), vbj = vr * (var.gstatcom()).vbinv(j);
    double fj = bj * bj / vbj;
    double pr = ftop(fj, 1, var.dfr());
    double ci = calmcci(mcfis, opt.cc(), 1, 1, var.dfr(), 1, vbj);
    if (j <= (var.gstatcom()).nx())
    { fout << "x" << left << setw(3) << j; }
    else
    {
      if ((var.gstatcom().tvec()).size() > 0) fout << left << setw(4) << "t";
      else fout << "t" << left << setw(3) << j - (var.gstatcom()).nx();
    }
    fout << right << setw(12) << bj;
    fout << setw(11) << sqrt(vbj) << setw(12) << bj - ci << setw(12) << bj + ci;
    fout << setw(12) << fj << setw(12) << pr << pind(pr) << "\n";
  }
  putline(fout, LENline);
}

void ancova2class::put_ancova2tab(ostream& fout,double sbc, int dfb, double sa, int dfa, double sbt, double saa, double sd, int dfd, double sr, int dfr, double st, int dft) const
// 共分散分析表を出力する
{
  fout << "\n";
  constexpr int LENline = 78;   // 区切り線の長さ
  string tstr;
  if (dfd > 0) tstr = "共分散分析表(ANCOVA2 table)";
  else tstr = "共分散分析表(非平行性プール)";
  int w = (LENline - strclen(tstr, CHARwide)) / 2 + tstr.size();
  fout << setw(w) << tstr;
  fout << "\n";
  fout << "要因              平方和    自由度    平均平方和         F値   有意確率p値\n";
  putline(fout, LENline);
  double v = 0.0, f = 0.0, pr = 0.0, fd = 0.0;
  double vr = calvfp(v, f, pr, sa, dfa, sr, dfr);
  tstr = "群差";
  int npad = tstr.size() - strclen(tstr, CHARwide);
  fout << left << setw(10 + npad) << tstr << right << setw(14) << sa;
  fout << setw(10) << dfa << setw(14) << v << setw(12) << f;
  fout << setw(14) << pr << pind(pr) << "\n";
  calvfp(v, f, pr, sbc, dfb, sr, dfr);
  tstr = "共通回帰";
  npad = tstr.size() - strclen(tstr, CHARwide);
  fout << left << setw(10 + npad) << tstr << right << setw(14) << sbc;
  fout << setw(10) << dfb << setw(14) << v << setw(12) << f;
  fout << setw(14) << pr << pind(pr) << "\n";
  putline(fout, LENline);
  calvfp(v, f, pr, saa, dfa, sr, dfr);
  tstr = "修正群差";
  npad = tstr.size() - strclen(tstr, CHARwide);
  fout << left << setw(10 + npad) << tstr << right << setw(14) << saa;
  fout << setw(10) << dfa << setw(14) << v << setw(12) << f;
  fout << setw(14) << pr << pind(pr) << "\n";
  calvfp(v, f, pr, sbt, dfb, sr, dfr);
  tstr = "全体回帰";
  npad = tstr.size() - strclen(tstr, CHARwide);
  fout << left << setw(10 + npad) << tstr << right << setw(14) << sbt;
  fout << setw(10) << dfb << setw(14) << v << setw(12) << f;
  fout << setw(14) << pr << pind(pr) << "\n";
  putline(fout, LENline);
  if (dfd > 0)
  {
    calvfp(v, f, pr, sd, dfd, sr, dfr);
    fd = f;
    tstr = "非平行性";
    npad = tstr.size() - strclen(tstr, CHARwide);
    fout << left << setw(10 + npad) << tstr << right << setw(14) << sd;
    fout << setw(10) << dfd << setw(14) << v << setw(12) << f;
    fout << setw(14) << pr << pind(pr) << "\n";
  }
  tstr = "残差";
  npad = tstr.size() - strclen(tstr, CHARwide);
  fout << left << setw(10 + npad) << tstr << right << setw(14) << sr;
  fout << setw(10) << dfr << setw(14) << vr << "\n";
  putline(fout, LENline);
  tstr = "全体";
  npad = tstr.size() - strclen(tstr, CHARwide);
  fout << left << setw(10 + npad) << tstr << right << setw(14) << st;
  fout << setw(10) << dft << endl;
}

void ancova2class::put_macmp(ostream& fout, double sr, int dfr, const varclass& var)
// 修正群差の多重比較を出力する
{
  constexpr int LENline = 77;   // 区切り線の長さ
  if (sr <= 0.0) return;        // 出力なし
  int g = (var.gsmap()).size(), d = (var.gstatcom()).dim() - 1;
  int df1 = 1, dc = 1;
  string method = "Tukey";
  switch (opt.ma())
  {
  case mcbon : method = "Bonferroni";  dc = g * (g - 1) / 2;
    break;
  case mcfis : method = "Fisher";
    break;
  case mcsch : method = "Scheffe";  df1 = g - 1;
  }
  fout << "\n・修正群差の" << method << "型多重比較:共変数=全体平均値  ";
  if (opt.ma() == mctuk)
  {
    fout << "群数=" << g << "  自由度=" << dfr << "\n";
    fout << "群 :  修正平均値 - 群 :  修正平均値    修正群差         q値   有意確率p値\n";
  }
  else
  {
    fout << "第1自由度=" << df1 << "  第2自由度=" << dfr << "\n";
    fout << "群 :  修正平均値 - 群 :  修正平均値    修正群差         F値   有意確率p値\n";
  }
  putline(fout, LENline);
  double vr = sr / dfr;
  vec<double> dxv(d - 1);
  double yc = 0.0;
  for (int k = 1; k < d; ++k) yc += (var.gstatcom()).b(k) * (var.gstatall()).me(k);
  mapsg_citr gitri = (var.gsmap()).begin();
  for (int i = 0; i < g - 1; ++i, ++gitri)
  {
    double ni = (double)(gitri->second.n());
    double bc0i = gitri->second.my();
    for (int k = 1; k < d; ++k) bc0i -= (var.gstatcom()).b(k) * gitri->second.me(k);
    double yi = bc0i + yc;
    mapsg_citr gitrj = gitri;
    ++gitrj;
    int j = i + 1;
    while (gitrj != (var.gsmap()).end())
    {
      double bc0j = gitrj->second.my();
      for (int k = 1; k < d; ++k)
      {
        bc0j -= (var.gstatcom()).b(k) * gitrj->second.me(k);
        dxv.assign(k - 1, gitri->second.me(k) - gitrj->second.me(k));
      }
      double yj = bc0j + yc;
      double dy = bc0j - bc0i;
      double vd = (1.0 / ni + 1.0 / (double)(gitrj->second.n()) + (var.gstatcom()).vxr(dxv)) * vr;
      double f = (sgn(vd) > 0) ? dy * dy / vd / (double)df1 : 0.0;
      double pr = calmcpr(opt.ma(), f, g, df1, dfr, dc);
      fout << left << setw(3) << j + 1 << ":" << right << setw(12) << yj << " - ";
      fout << left << setw(3) << i + 1 << ":" << right << setw(12) << yi << setw(12) << dy;
      fout << setw(12) << f << setw(14) << pr << pind(pr) << "\n";
      ++gitrj;  ++j;
    }
  }
  putline(fout, LENline);
}

void ancova2class::put_maci(ostream& fout, double sr, int dfr, const varclass& var)
// 修正群差の信頼区間を出力する
{
  const int LENline = 57;       // 区切り線の長さ
  if (sr <= 0.0) return;        // 出力なし
  int g = (var.gsmap()).size(), d = (var.gstatcom()).dim() - 1;
  int df1 = 1, dc = 1;
  string method = "Tukey";
  switch (opt.ma())
  {
  case mcbon : method = "Bonferroni";  dc = g * (g - 1) / 2;
    break;
  case mcfis : method = "Fisher";
    break;
  case mcsch : method = "Scheffe";  df1 = g - 1;
  }
  fout << "\n・修正群差の" << method << "型" << opt.cc() * 100.0 << "%同時信頼区間(simultaneous confidence interval)\n";
  fout << "群  - 群     修正群差    絶対精度        下限        上限\n";
  putline(fout, LENline);
  double vr = sr / dfr;
  vec<double> dxv(d - 1);
  mapsg_citr gitri = (var.gsmap()).begin();
  for (int i = 0; i < g - 1; ++i, ++gitri)
  {
    double ni = (double)(gitri->second.n());
    double bc0i = gitri->second.my();
    for (int k = 1; k < d; ++k) bc0i -= (var.gstatcom()).b(k) * gitri->second.me(k);
    mapsg_citr gitrj = gitri;
    ++gitrj;
    int j = i + 1;
    while (gitrj != (var.gsmap()).end())
    {
      double bc0j = gitrj->second.my();
      for (int k = 1; k < d; ++k)
      {
        bc0j -= (var.gstatcom()).b(k) * gitrj->second.me(k);
        dxv.assign(k - 1, gitri->second.me(k) - gitrj->second.me(k));
      }
      double dy = bc0j - bc0i;
      double vd = (1.0 / ni + 1.0 / gitrj->second.n() + (var.gstatcom()).vxr(dxv)) * vr;
      double ci = calmcci(opt.ma(), opt.cc(), g, df1, dfr, dc, vd);
      fout << left << setw(4) << j + 1 << "- " << setw(3) << i + 1 << right;
      fout << setw(12) << dy << setw(12) << ci << setw(12) << dy - ci << setw(12) << dy + ci << "\n";
      ++gitrj;  ++j;
    }
  }
  putline(fout, LENline);
}

void ancova2class::put_mdcmp(ostream& fout, double sr, int dfr, const varclass& var)
// 非平行性の多重比較を出力する
{
  const int LENline = 57;       // 区切り線の長さ
  if (sr <= 0.0) return;        // 出力なし
  int g = (var.gsmap()).size(), d = (var.gstatcom()).dim() - 1;
  int df1 = d - 1, dc = dc = g * (g - 1) / 2;
  string method = "Bonferroni";
  switch (opt.md())
  {
  case mcfis : method = "Fisher";  dc = 1;
    break;
  case mcsch : method = "Scheffe";  df1 = (g - 1) * df1;  dc = 1;
  }
  fout << "\n・非平行性の" << method << "型多重比較:第1自由度=" << df1 << "  第2自由度=" << dfr << "\n";
  fout << "群  - 群   回帰係数差平方和         F値   有意確率p値\n";
  putline(fout, LENline);
  double vr = sr / dfr;
  vec<double> dbv(d - 1);
  mapsg_citr gitri = (var.gsmap()).begin();
  for (int i = 0; i < g - 1; ++i, ++gitri)
  {
    mapsg_citr gitrj = gitri;
    ++gitrj;
    int j = i + 1;
    while (gitrj != (var.gsmap()).end())
    {
      gstat gstatc(opt, (var.gstatcom()).p());
      gstatc += gitri->second;  gstatc += gitrj->second;
      gstatc.calreg();
      for (int k = 1; k < d; ++k)
      { dbv.assign(k - 1, gitrj->second.b(k) - gitri->second.b(k)); }
      double sdb = gstatc.vxr(dbv);
      double f = (sgn(vr) > 0) ? sdb / vr / df1 : 0.0;
      double pr = calmcpr(opt.md(), f, g, df1, dfr, dc);
      fout << left << setw(4) << j + 1 << "- " << setw(3) << i + 1 << right;
      fout << setw(18) << sdb << setw(12) << f;
      fout << setw(14) << pr << pind(pr) << "\n";
      ++gitrj;  ++j;
    }
  }
  putline(fout, LENline);
}

int ancova2class::get(const char* c_fname) // データを読み込み保存・累積する
{
  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 ancova2class::put(int nc)  // 集計結果を出力する
{
  if (opt.res("res") == false) return;
  varclass var(opt, id.ni());     // 変数オブジェクト
  var.cal(opt, datalst);
  if (var.nt() < 1)             // 有効なデータなし
  {
    put_resid(cout);
    cout << "※有効なデータなし！\n\n" << endl;
    return;
  }
  double vd = 0.0, fd = 0.0, pr = 0.0;
  calvfp(vd, fd, pr, var.sd(), var.dfd(), var.sr(), var.dfr());
  put_resid(cout);
  put_mreg(cout, var);
  if (opt.res("bt")) put_bt(cout, var);
  if (var.ng() < 2)
  { cout << "※共分散分析不可能！\n"; }
  else
  {
    if (opt.res("atab")) put_ancova2tab(cout, var.sbc(), var.dfb(), var.sa(), var.dfa(), var.sbt(), var.saa(), var.sd(), var.dfd(), var.sr(), var.dfr(), var.st(), var.dft());
    if (opt.res("ma"))
    {
      put_macmp(cout, var.sr(), var.dfr(), var);
      put_maci(cout, var.sr(), var.dfr(), var);
    }
    if (opt.res("md")) put_mdcmp(cout, var.sr(), var.dfr(), var);
    if ((opt.res("def") && fd < 1.0) || opt.res("pd"))
    {
      double srt = var.st() - var.saa() - var.sbt();
      int dfrt = var.dft() - var.dfa() - var.dfb();
      put_ancova2tab(cout, var.sbc(), var.dfb(), var.sa(), var.dfa(), var.sbt(), var.saa(), 0.0, 0, srt, dfrt, var.st(), var.dft());
      if (opt.res("ma"))
      {
        put_macmp(cout, srt, dfrt, var);
        put_maci(cout, srt, dfrt, var);
      }
    }
  }
  cout << "\n" << endl;
}

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