// crosov.C Version 7.44  Copyleft (c) 杉本典夫
// Time-stamp: <2025-10-28 17:31:20 tac>
//
// クロスオーバー・デザインによる分散分析。
// 主として生物学的同等性試験解析用。
//
// 使用ライブラリ:libdans.so, libfun.so
// コンパイル法 :g++ -O -s -W crosov crosov.C
//             -I/usr/local/include -L/usr/local/lib -ldans -lfun

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

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

const vecstr USAGE =            // プログラム使用法表示用
{
  "crosov  Version 7.44         Copyleft(Vo\\oV) 2025 杉本典夫",
  "書式:crosov [<オプション> …]  [<データファイル> …]",
  "機能:最初の項目を群項目、以後の項目を時期項目としてクロスオーバーデザイン",
  "     による分散分析をするのだ。群数と時期数は同じでなければいかんのだ。",
  "<オプション>",
  "  -al[=<有意水準α>]:有意水準αを指定するのだ。(0＜α＜1,-al:-al=0.05)",
  "  -cc[=<信頼係数>]:信頼係数を指定するのだ。(0≦信頼係数<1)",
  "    -cc:-cc=0.95(オプション省略時)",
  "  -dl[=<検出差>]:検出差の指定。(対数変換後の実測値,-dl:-dl=log(1.25))",
  "  -dp[=<検出差>]:検出差の指定。(薬剤1の平均値に対する比率,-dp:-dp=0.2)",
  "  -dg=<処置名1>,<処置名2>[,<処置名3>,…]:処置名を指定するのだ。",
  "  -f=<ファイル名>:指定されたファイルからオプションを読み込むのだ。",
  "  -g=<群コード1>,<群コード2>[,<群コード3>,…]:群のコードを指定するのだ。",
  "  -h[=<コマンド名>]:オンラインマニュアルを表示するのだ。(-h:-h=man)",
  "  -log={10|e}:対数の種類を指定するのだ。(e:自然対数(省略時) 10:常用対数)",
  "  -po[=<検出力(1-β)>]:検出力(1-β)を指定するのだ。(0＜1-β＜1,-po:-po=0.8)",
  "  -res=[-][{def|all|at|be|pa|stat|tt|tt2|ttw|ut|utc|uto}]",
  "      [,…]:解析結果の出力項目を指定するのだ。",
  "    -:-を付けた項目を出力しない(-res=-:全ての項目を出力しない)",
  "    def:オプション省略時と同じ項目を出力する  all:全ての項目を出力する",
  "    at:分散分析  be:生物的同等性試験用設定  pa:検出力分析  stat:基礎統計量",
  "    tt:自動判定t検定  tt2:2標本t検定  ttw:Welchの検定  ut:自動判定U検定",
  "    utc:連続修正有U検定  uto:連続修正無U検定",
  "    オプション省略時:データに合わせて自動判定"
};

const string Title = "=== クロスオーバー用分散分析 ===";
constexpr double DLdef = 1.25;  // 省略時検出差(対数変換前の値)
constexpr double DPdef = 0.2;   // 省略時検出差(％)

using vecstatv = vector<statv<double>>; // 基礎統計量ベクトルの短名定義

enum sumtype                    // 度数集計の種類
{
  sumnon,                       // 集計なし
  sumorder,                     // 順序効果用集計
  sumperiod,                    // 時期効果用集計
  sumtreat,                     // 処置効果用集計
  sumetc                        // その他
};

enum logtype                    // 対数の種類
{
  logtnon,                      // 指定なし
  logte,                        // 自然対数
  logt10,                       // 常用対数
  logtetc                       // その他
};

class option                    // オプションクラス
{
private:
  string v_path;                // プログラムのパス
  string v_pname;               // プログラム名
  double v_al;                  // 有意水準
  double v_cc;                  // 信頼係数
  double v_dl;                  // 検出差
  double v_dp;                  // 検出差(比率)
  logtype v_logt;               // 対数の種類フラッグ
  double v_po;                  // 検出力
  vecstr dgvec;                 // 処置名ベクトル
  vecstr gvec;                  // 群コードベクトル
  lstr reslst;                  // 分析結果出力リスト
  mapsb resmap;                 // 分析結果出力の有無マップ
  testtype v_tt;                // t検定の種類
  testtype v_ut;                // U検定の種類
public:
  option(const char* c_pname = nullptr); // コンストラクタ
  ~option()                              // デストラクタ
  {}
  const string& path() const    // プログラムのパスを返す
  { return v_path; }
  const string& pname() const   // プログラム名を返す
  { return v_pname; }
  double al() const             // 有意水準を返す
  { return v_al; }
  double cc() const             // 信頼係数を返す
  { return v_cc; }
  double dl() const             // 検出差を返す
  { return v_dl; }
  double dp() const             // 検出差(%)を返す
  { return v_dp; }
  logtype logt() const          // 対数の種類を返す
  { return v_logt; }
  double po() const             // 検出力を返す
  { return v_po; }
  const string& dg(int idx) const // 処置名を返す
  { return dgvec[idx]; }
  const string& g(int idx) const // 群コードを返す
  { return gvec[idx]; }
  bool res(const string& code); // 分析結果出力フラッグを返す
  testtype tt() const           // t検定の種類を返す
  { return v_tt; }
  testtype ut() const           // U検定の種類を返す
  { return v_ut; }
  void setdef()                 // 省略時オプション指定を設定する
  { getoptdef(*this); }
  void set(const string& optstr, char devc = OPTdev);
  // オプション指定によってオプションを設定する
  int check(int g);             // オプションのチェックとセット
  int set_res();                // -resオプションのセット
  int search_g(const string& gcode) const; // 群コードを検索する
};

option::option(const char* c_pname) : v_path(), v_pname(), v_al(ALdef), v_cc(CCdef), v_dl(DLdef), v_dp(DPdef), v_logt(logte), v_po(POdef), dgvec(), gvec(), reslst(), resmap(), v_tt(testaut), v_ut(testaut)
// コンストラクタ
{
  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 respos = resmap.find(code);
  if (respos == 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] == "al")           // 有意水準指定
  {
    v_al = opt[1].empty() ? ALdef : strtodf(opt[1]);
    if (v_al <= 0.0 || 1.0 <= v_al) err = badopt(pname(), optstr, USAGE, 0.0, 1.0);
  }
  else 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] == "dl")      // 検出差の指定
  {
    if (opt[1].empty()) v_dl = (v_logt == logt10) ? log10(DLdef) : log(DLdef);
    else v_dl = strtodf(opt[1]);
    if (v_dl < 0.0 || FDmax < v_dl) err = badopt(pname(), optstr, USAGE, 0.0, FDmax);
  }
  else if (opt[0] == "dp")      // 検出差(%)の指定
  {
    v_dp = opt[1].empty() ? DPdef : strtodf(opt[1]);
    if (v_dp < 0.0 || FDmax < v_dp) err = badopt(pname(), optstr, USAGE, 0.0, FDmax);
  }
  else if (opt[0] == "dg")      // 処置名の指定
  {
    if (strsplit(dgvec, opt[1], CHRdelm) < 1) err = badopt(pname(), optstr, USAGE);
  }
  else if (opt[0] == "f")       // オプション指定ファイル
  { err = fgetopt(*this, opt[1], true); }
  else if (opt[0] == "g")       // 群コードの指定
  {
    if (strsplit(gvec, opt[1], CHRdelm) < 1) err = badopt(pname(), optstr, USAGE);
  }
  else if (opt[0] == "h")       // オンラインヘルプ表示
  { err = help(pname(), USAGE, opt[1]); }
  else if (opt[0] == "log")     // 対数の種類の指定
  {
    if (opt[1].empty() || opt[1] == "e") v_logt = logte;
    else if (opt[1] == "10") v_logt = logt10;
    else err = badopt(pname(), optstr, USAGE);
  }
  else if (opt[0] == "po")      // 検出力の指定
  {
    v_po = opt[1].empty() ? POdef : strtodf(opt[1]);
    if (v_po <= 0.0 || 1.0 <= v_po) err = badopt(pname(), optstr, USAGE, 0.0, 1.0);
  }
  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(int g)        // オプションのチェック
{
  if (g != 2) return disperr(pname(), "群数は2つだけなのだ！");
  if ((int)(dgvec.size()) < g)  // 処置名が足らない場合
  {
    int i = dgvec.size();
    dgvec.resize(g);
    while (i < g)
    {
      ostringstream ostr;
      ostr << "DG" << i + 1;
      dgvec[i] = ostr.str();
      ++i;
    }
  }
  if (g < (int)(dgvec.size())) return disperr(pname(), "集計項目数が不適なのだ！");
  if ((int)(gvec.size()) < g)   // 群コードが足らない場合
  {
    int i = gvec.size();
    gvec.resize(g);
    while (i < g)
    {
      ostringstream ostr;
      ostr << i + 1;
      gvec[i] = ostr.str();
      ++i;
    }
  }
  if (set_res() < 0) return -2; // -resオプションのセット
  return 0;                     // エラーなし
}

int option::set_res()           // -resオプションのセット
{
  resmap["res"] = resmap["def"] = resmap["at"] = resmap["stat"] = true;
  resmap["be"] = resmap["pa"] = resmap["tt"] = resmap["ut"] = false;
  if (reslst.size() > 0)
  {
    resmap["def"] = resmap["at"] = resmap["stat"] = false;
    for (const string& resstr : reslst)
    {
      mapsb_citr respos = resmap.find(resstr);
      if (resstr == "-")        // 全て出力せず
      { resmap["res"] = false; }
      else if (resstr[0] == '-' && resstr.size() > 1) // 特定の項目を出力せず
      {
        respos = resmap.find(resstr.substr(1));
        if (respos == resmap.end())
        { return disperr(pname(), "-res=" + resstr + "オプション指定が不適当！"); }
        else
        { resmap[resstr.substr(1)] = false; }
      }
      else if (resstr == "def") // 省略時出力
      { resmap["def"] = resmap["at"] = resmap["stat"] = true; }
      else if (resstr == "all")
      {
        resmap["at"] = resmap["be"] = resmap["pa"] = resmap["stat"] = resmap["tt"] = resmap["ut"] = true;
        v_al = 0.1;  v_cc = 0.9;  v_dp = 0.2;  v_po = 0.8;
        v_dl = (v_logt == logt10) ? log10(DLdef) : log(DLdef);
      }
      else if (resstr == "be")  // 生物学的同等性試験用設定
      {
        v_al = 0.1;  v_cc = 0.9;  v_dp = 0.2;  v_po = 0.8;
        v_dl = (v_logt == logt10) ? log10(DLdef) : log(DLdef);
        resmap["pa"] = true;
      }
      else if (resstr == "tt2")
      { resmap["tt"] = true;  v_tt = testorg; } // 2標本t検定
      else if (resstr == "ttw")
      { resmap["tt"] = true;  v_tt = testcor; } // Welchの検定
      else if (resstr == "utc")
      { resmap["ut"] = true;  v_ut = testcor; } // 連続修正有り
      else if (resstr == "uto")
      { resmap["ut"] = true;  v_ut = testorg; } // 連続修正無し
      else if (respos == resmap.end())
      { return disperr(pname(), "-res=" + resstr + "オプション指定が不適当！"); }
      else
      { resmap[resstr] = true; }
    }
  }
  return 0;                     // エラーなし
}

int option::search_g(const string& gcode) const // 群コードを検索する
{
  bool hit = false;
  int g = gvec.size(), idx = 0;
  while (hit == false && idx < g)
  {
    if (gcode == gvec[idx]) hit = true;
    else ++idx;
  }
  if (hit) return idx;          // 発見
  return -1;
}

class varclass                  // 変数クラス
{
private:
  vector<vecstatv> gsvec;       // 群別統計量2重ベクトル
  set<double> v_sumset;         // クロス前後の和データセット
  set<double> v_difset;         // クロス前後の差データセット
  vecstatv tgstatvec;           // 群統計量ベクトル
  vecstatv tpstatvec;           // 時期統計量ベクトル
  vecstatv tdstatvec;           // 薬剤統計量ベクトル
  statv<double> tstat;          // 全体統計量
  double v_sg;                  // 群(順序)平方和
  int v_dfg;                    // 群(順序)自由度
  double v_sbr;                 // 群(順序)残差平方和
  int v_dfbr;                   // 群(順序)残差自由度
  double v_ss;                  // 個体平方和
  int v_dfs;                    // 個体自由度
  double v_sp;                  // 時期平方和
  double v_sd;                  // 処置平方和
  double v_sr;                  // 残差平方和
  int v_dfr;                    // 残差自由度
  int v_dft;                    // 全体自由度
  restype v_resvar;             // 処理結果の種類
  void cal_crosov();        // クロスオーバー用各種パラメータを計算する
public:
  varclass(size_t g0) : gsvec(g0), v_sumset(), v_difset(), tgstatvec(g0), tpstatvec(g0), tdstatvec(g0), tstat(), v_sg(0.0), v_dfg(0), v_sbr(0.0), v_dfbr(0), v_ss(0.0), v_dfs(0), v_sp(0.0), v_sd(0.0), v_sr(0.0), v_dfr(0), v_dft(0), v_resvar(resnon)
  // デフォールトコンストラクタ
  {}
  ~varclass()                   // デストラクタ
  {}
  size_t g() const              // 群数を返す
  { return gsvec.size(); }
  int tgn(size_t idx) const     // 群別例数を返す
  { return (idx < tgstatvec.size()) ? tgstatvec[idx].n() : 0; }
  double tgmean(size_t idx) const // 群別平均値を返す
  { return (idx < tgstatvec.size()) ? tgstatvec[idx].mean() : 0.0; }
  double tgvr(size_t idx) const // 群別不偏分散を返す
  { return (idx < tgstatvec.size()) ? tgstatvec[idx].vr() : 0.0; }
  int tpn(size_t idx) const     // 時期別例数を返す
  { return (idx < tpstatvec.size()) ? tpstatvec[idx].n() : 0; }
  double tpmean(size_t idx) const // 時期別平均値を返す
  { return (idx < tpstatvec.size()) ? tpstatvec[idx].mean() : 0.0; }
  double tpvr(size_t idx) const // 時期別不偏分散を返す
  { return (idx < tpstatvec.size()) ? tpstatvec[idx].vr() : 0.0; }
  int tdn(size_t idx) const     // 処置別例数を返す
  { return (idx < tdstatvec.size()) ? tdstatvec[idx].n() : 0; }
  double tdmean(size_t idx) const // 処置別平均値を返す
  { return (idx < tdstatvec.size()) ? tdstatvec[idx].mean() : 0.0; }
  double tdvr(size_t idx) const // 処置別不偏分散を返す
  { return (idx < tdstatvec.size()) ? tdstatvec[idx].vr() : 0.0; }
  int n(size_t iidx, size_t jidx) const // 群別時期別例数を返す
  { return (iidx < gsvec.size() && jidx < gsvec.size()) ? (gsvec[iidx])[jidx].n() : 0; }
  double mean(size_t iidx, size_t jidx) const // 群別時期別平均値を返す
  { return (iidx < gsvec.size() && jidx < gsvec.size()) ? (gsvec[iidx])[jidx].mean() : 0.0; }
  double vr(size_t iidx, size_t jidx) const // 群別時期別不偏分散を返す
  { return (iidx < gsvec.size() && jidx < gsvec.size()) ? (gsvec[iidx])[jidx].vr() : 0.0; }
  int sn(size_t idx) const      // 群別クロス和例数を返す
  { return (idx < gsvec.size()) ? (gsvec[idx])[2].n() : 0; }
  double smean(size_t idx) const // 群別クロス和平均値を返す
  { return (idx < gsvec.size()) ? (gsvec[idx])[2].mean() : 0.0; }
  double sss(size_t idx) const // 群別クロス和平方和を返す
  { return (idx < gsvec.size()) ? (gsvec[idx])[2].ss() : 0.0; }
  double svr(size_t idx) const // 群別クロス和不偏分散を返す
  { return (idx < gsvec.size()) ? (gsvec[idx])[2].vr() : 0.0; }
  int dn(size_t idx) const      // 群別クロス差例数を返す
  { return (idx < gsvec.size()) ? (gsvec[idx])[3].n() : 0; }
  double dmean(size_t idx) const // 群別クロス差平均値を返す
  { return (idx < gsvec.size()) ? (gsvec[idx])[3].mean() : 0.0; }
  double dss(size_t idx) const // 群別クロス差平方和を返す
  { return (idx < gsvec.size()) ? (gsvec[idx])[3].ss() : 0.0; }
  double dvr(size_t idx) const // 群別クロス差不偏分散を返す
  { return (idx < gsvec.size()) ? (gsvec[idx])[3].vr() : 0.0; }
  int tn() const                // 全体の例数を返す
  { return tstat.n(); }
  double tmean() const          // 全体の平均値を返す
  { return tstat.mean(); }
  double tvr() const            // 全体の不偏分散を返す
  { return tstat.vr(); }
  double sg() const             // 群(順序)平方和を返す
  { return v_sg; }
  int dfg() const               // 群(順序)自由度を返す
  { return v_dfg; }
  double sbr() const            // 群(順序)残差平方和を返す
  { return v_sbr; }
  int dfbr() const              // 群(順序)残差自由度を返す
  { return v_dfbr; }
  double ss() const             // 個体平方和を返す
  { return v_ss; }
  int dfs() const               // 個体自由度を返す
  { return v_dfs; }
  double sp() const             // 時期平方和を返す
  { return v_sp; }
  double sd() const             // 処置平方和を返す
  { return v_sd; }
  double sr() const             // 残差平方和を返す
  { return v_sr; }
  int dfr() const               // 残差自由度を返す
  { return v_dfr; }
  double vr() const             // 残差分散を返す
  { return (dfr() > 0) ? sr() / (double)(dfr()) : 0.0; }
  int dft() const               // 全体自由度を返す
  { return v_dft; }
  double st() const             // 全体の平方和を返す
  { return tstat.ss(); }
  restype resvar() const        // 処理結果の種類を返す
  { return v_resvar; }
  const set<double>& sumset() const // クロス前後の和データセットを返す
  { return v_sumset; }
  const set<double>& difset() const // クロス前後の差データセットを返す
  { return v_difset; }
  void cal(option& opt, ldata& datalst); // データを集計する
};

void varclass::cal_crosov() // クロスオーバー用各種パラメータを計算する
{
  size_t g = gsvec.size();
  v_resvar = resok;
  for (size_t i = 0; i < g; ++i)
  {
    for (size_t j = 0; j < g; ++j)
    {
      tgstatvec[i] += (gsvec[i])[j]; // 群データ累積
      tpstatvec[j] += (gsvec[i])[j]; // 時期データ累積
      size_t k = i + j;
      if (k >= g) k -= g;
      tdstatvec[k] += (gsvec[i])[j]; // 薬剤データ累積
      if (gsvec[i][j].n() < 1) v_resvar = reserr; // 分散分析不可能
    }
  }
  statv<double> gstat;
  statv<double> dstat;
  for (size_t i = 0; i < g; ++i)
  {
    v_sg += (gsvec[i])[2].mean() * (gsvec[i])[2].mean() * (gsvec[i])[2].n();
    v_sd += (gsvec[i])[3].mean() * (gsvec[i])[3].mean() * (gsvec[i])[3].n();
    gstat += (gsvec[i])[2];  dstat += (gsvec[i])[3];
    tstat += tgstatvec[i];
  }
  double cf = tstat.mean() * tstat.mean() * tstat.n();
  v_ss -= cf;
  cf = gstat.mean() * gstat.mean() * gstat.n();
  v_sg = (v_sg  - cf) * 0.5;  v_sp = v_sd;
  cf = dstat.mean() * dstat.mean() * dstat.n();
  v_sd = (v_sd - cf) * 0.5;
  double mt = ((gsvec[0])[3].sum() - (gsvec[1])[3].sum()) / dstat.n();
  cf = mt * mt * dstat.n();
  v_sp = (v_sp - cf) * 0.5;
  v_sbr = v_ss - v_sg;
  v_sr = ((gsvec[0])[3].ss() + (gsvec[1])[3].ss()) * 0.5;
  if (v_ss < FDeps) v_ss = 0.0;
  if (v_sg < FDeps) v_sg = 0.0;
  if (v_sp < FDeps) v_sp = 0.0;
  if (v_sd < FDeps) v_sd = 0.0;
  if (v_sr < FDeps) v_sr = 0.0;
  v_dft = tstat.n() - 1;  v_dfg = g - 1;  v_dfs = gstat.n() - 1;
  v_dfbr = v_dfs - v_dfg;  v_dfr = v_dft - v_dfs - v_dfg * 2;
}

void varclass::cal(option& opt, ldata& datalst)
// データを集計する
{
  size_t g = gsvec.size();
  for (size_t i = 0; i < g; ++i) gsvec[i].resize(g + 2);
  vecdbl fdatavec(g);
  int nc = 0;
  v_ss = 0.0;                   // 個体平方和
  for (auto& cdatai : datalst)
  {
    int idx = opt.search_g(cdatai.sdata(0));
    if (idx >= 0)               // 群コード該当
    {
      bool blank = false;
      size_t j = 0;
      while (blank == false && j < g)
      {
        fdatavec[j] = cdatai.fdata(j + 1);
        if (fdatavec[j] == BLANKf) blank = true;
        else ++j;
      }
      if (blank == false)
      {
        double sum = fdatavec[0] + fdatavec[1]; // クロス前後の和
        double dif = fdatavec[1] - fdatavec[0]; // クロス前後の差
        (gsvec[idx])[0].sigma(fdatavec[0]);
        (gsvec[idx])[1].sigma(fdatavec[1]);
        (gsvec[idx])[2].sigma(sum);
        (gsvec[idx])[3].sigma(-dif); // 符号を逆転した前後差を累積する
        v_sumset.insert(sum);        // 新規データだけ挿入する
        v_difset.insert(dif);
        if (idx == 1) v_difset.insert(-dif); // 群2は符合を逆転した差も挿入する
        v_ss += sum * sum / 2.0; // 個体平方和用の値を累積する
        ++nc;
      }
    }
  }
  cal_crosov();
}

class crosovclass               // crosovメインクラス
{
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_coatab(ostream& fout, const varclass& var) const;
  // クロスオーバー用分散分析表を出力する
  void put_dgest(ostream& fout, option& opt, const varclass& var);
  // 処置に関する検定特性を出力する
  void put_orderest(ostream& fout, option& opt, const varclass& var) const;
  // 順序効果の信頼区間を出力する
  void put_dgd(ostream& fout, option& opt, int iidx, int ni, double mdi, double dij, double sed, double ci, double sedm, double dmin) const;
  // 処置差とその信頼区間と検出差を出力する
  void put_logdn(ostream& fout, option& opt, double ta, double za, double zb, int ni, double dij, double sed, int dfr, double vr) const;
  // 対数検出差に関する検出力と必要例数を出力する
  void put_prodn(ostream& fout, option& opt, double ta, double za, double zb, int ni, double mdi, double dij, double sed, int dfr, double vr) const;
  // 割合検出差に関する検出力と必要例数を出力する
  int cal_nci(option& opt, double& ci, ldbl& po0, double dd, double d, double vr) const;
  // 信頼区間による必要例数を計算して返す
  void put_ttest(ostream& fout, option& opt, const varclass& var) const;
  // クロスオーバー用t検定結果を出力する
  void put_tstat(ostream& fout, option& opt, string& tstr, int n1, double m1, double ss1, double v1, int n2, double m2, double ss2, double v2) const;
  // 2標本t検定の基礎統計量と検定結果を出力する
  void put_mdci(ostream& fout, option& opt, double md, int df, double vmd) const;
  // 平均値の差の信頼区間を出力する
  void put_utest(ostream& fout, option& opt, const varclass& var);
  // クロスオーバー用U検定結果を出力する
  void put_umrdci(ostream& fout, option& opt, ftable& ftab, double u, double vu) const;
  // 順位平均値の差の信頼区間とU値の信頼区間を出力する
  void sum_ftab(option& opt, ftable& ftab, sumtype stype, const varclass& var);
  // クロスオーバー用U検定用用にクロス前後の和／差／クロス差を度数集計する
  int search_data(const set<double>& targset, double fdata) const;
  // 数字データをサーチして番号を返す
public:
  crosovclass() : fname(), id(), datalst() // コンストラクタ
  {}
  ~crosovclass()                // デストラクタ
  {}
  int get(const char* c_fname, option& opt); // データを読み込み保存・累積する
  void put(option& opt);                     // 集計結果を出力する
};

void crosovclass::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 << "時期項目" << setw(2) << i;
    fout << ":" << id.item(i) << "\n";
  }
}

void crosovclass::put_stat(ostream& fout, option& opt, const varclass& var) const
// 基礎統計量を出力する
{
  constexpr int LENline = 77;   // 区切り線の長さ
  constexpr int LENhed = 9;     // ヘッダの長さ
  size_t g = var.g();
  fout << "\n群       :群別・全時期合計基礎統計量\n";
  putline(fout, LENline);
  for (size_t i = 0; i < g; ++i)
  {
    fout << left << setw(LENhed) << opt.g(i) << right << ":";
    putmsd(fout, var.tgn(i), var.tgmean(i), var.tgvr(i));
  }
  putline(fout, LENline);
  fout << "\n時期     :時期別・全群合計基礎統計量\n";
  putline(fout, LENline);
  for (size_t i = 0; i < g; ++i)
  {
    fout << left << setw(9) << i + 1 << right << ":";
    putmsd(fout, var.tpn(i), var.tpmean(i), var.tpvr(i));
  }
  putline(fout, LENline);
  fout << "\n処置     :処置別・全時期合計基礎統計量\n";
  putline(fout, LENline);
  for (size_t i = 0; i < g; ++i)
  {
    int npad = opt.dg(i).size() - strclen(opt.dg(i), CHARwide);
    fout << left << setw(LENhed + npad) << opt.dg(i) << right << ":";
    putmsd(fout, var.tdn(i), var.tdmean(i), var.tdvr(i));
  }
  putline(fout, LENline);
  fout << "\n群   時期:群別・時期別基礎統計量\n";
  putline(fout, LENline);
  for (size_t i = 0; i < g; ++i)
  {
    for (size_t j = 0; j < g; ++j)
    {
      int npad = opt.g(i).size() - strclen(opt.g(i), CHARwide);
      fout << left << setw(LENhed - 4 + npad) << opt.g(i) << right;
      fout << setw(4) << j + 1 << ":";
      putmsd(fout, var.n(i, j), var.mean(i, j), var.vr(i, j));
    }
  }
  putline(fout, LENline);
  fout << "全体     :";
  putmsd(fout, var.tn(), var.tmean(), var.tvr());
}

void crosovclass::put_coatab(ostream& fout, const varclass& var) const
// クロスオーバー用分散分析表を出力する
{
  constexpr int LENline = 78;   // 区切り線の長さ
  fout << "\n";
  string tstr("分散分析表(ANOVA table)");
  int npad = tstr.size() - strclen(tstr, CHARwide);
  fout << setw(51 + npad) << tstr << "\n";
  fout << "要因              平方和    自由度    平均平方和         F値   有意確率p値\n";
  putline(fout, LENline);
  double v = 0.0, f = 0.0, pr = 0.0;
  double vr = calvfp(v, f, pr, var.sg(), var.dfg(), var.sbr(), var.dfbr());
  tstr = "順序";
  npad = tstr.size() - strclen(tstr, CHARwide);
  fout << left << setw(10 + npad) << tstr << right << setw(14) << var.sg();
  fout << setw(10) << var.dfg() << 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) << var.sbr();
  fout << setw(10) << var.dfbr() << setw(14) << vr << "\n";
  putline(fout, LENline);
  calvfp(v, f, pr, var.ss(), var.dfs(), var.sr(), var.dfr());
  tstr = "観察個体";
  npad = tstr.size() - strclen(tstr, CHARwide);
  fout << left << setw(10 + npad) << tstr << right << setw(14) << var.ss();
  fout << setw(10) << var.dfs() << setw(14) << v << setw(12) << f;
  fout << setw(14) << pr << pind(pr) << "\n";
  calvfp(v, f, pr, var.sp(), var.dfg(), var.sr(), var.dfr());
  tstr = "時期";
  npad = tstr.size() - strclen(tstr, CHARwide);
  fout << left << setw(10 + npad) << tstr << right << setw(14) << var.sp();
  fout << setw(10) << var.dfg() << setw(14) << v << setw(12) << f;
  fout << setw(14) << pr << pind(pr) << "\n";
  vr = calvfp(v, f, pr, var.sd(), var.dfg(), var.sr(), var.dfr());
  tstr = "処置";
  npad = tstr.size() - strclen(tstr, CHARwide);
  fout << left << setw(10 + npad) << tstr << right << setw(14) << var.sd();
  fout << setw(10) << var.dfg() << 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) << var.sr();
  fout << setw(10) << var.dfr() << setw(14) << vr << "\n";
  putline(fout, LENline);
  tstr = "全体";
  npad = tstr.size() - strclen(tstr, CHARwide);
  fout << left << setw(10 + npad) << tstr << right << setw(14) << var.st();
  fout << setw(10) << var.dft() << endl;
}

void crosovclass::put_dgest(ostream& fout, option& opt, const varclass& var)
// 処置に関する検定特性を出力する
{
  constexpr int LENline = 78;   // 区切り線の長さ
  size_t g = var.g();
  int nt = var.tdn(0);
  double za = ptoz(opt.al()), zb = ptoz((1.0 - opt.po()) * 2.0);
  double ta = ptot(opt.al(), var.dfr()), tb = ptot((1.0 - opt.po()) * 2.0);
  double tc = ptot(1.0 - opt.cc(), var.dfr());
  fout << "\n○検定特性:有意水準(α)=" << opt.al() * 100.0 << "%  検出力(1-β)=" << opt.po() * 100.0 << "%  信頼係数=" << opt.cc() * 100.0 << "%\n";
  putline(fout, LENline);
  put_orderest(fout, opt, var);
  for (size_t i = 0; i < g - 1; ++i)
  {
    int ni = var.n(i, 0);
    double mdi = var.tdmean(i);
    double mddi = var.dmean(i);
    for (size_t j = i + 1; j < g; ++j)
    {
      fout << "・処置効果：" << opt.dg(i) << " 対 " << opt.dg(j);
      fout << " (差=" << opt.dg(j) << " - " << opt.dg(i) << ")\n";
      int nj = var.n(j, 0);
      double mdj = var.tdmean(j);
      double mddj = var.dmean(j);
      double dij = (mddj - mddi) * 0.5;
      double sed = sqrt(var.vr() / (double)g * (1.0 / (double)ni + 1.0 / (double)nj));
      double ci = tc * sed;
      double sedm = (ni > 1 && nj > 1) ? sqrt(var.vr() / (double)g * (1.0 / (double)(ni - 1) + 1.0 / (double)(nj - 1))) : 0.0;
      double dmin = (ta + tb) * sedm;
      put_dgd(fout, opt, i, ni, mdi, dij, sed, ci, sedm, dmin);
      if (opt.res("pa"))
      {
        if (opt.dl() > 0.0) put_logdn(fout, opt, ta, za, zb, ni, dij, sed, var.dfr(), var.vr());
        if (opt.dp() > 0.0) put_prodn(fout, opt, ta, za, zb, ni, mdi, dij, sed, var.dfr(), var.vr());
      }
      putline(fout, LENline);
    }
  }
}

void crosovclass::put_orderest(ostream& fout, option& opt, const varclass& var) const
// 順序効果の信頼区間を出力する
{
  int n1 = var.sn(0), n2 = var.sn(1);
  int nt = n1 + n2, df = n1 + n2 - 2;
  double ss1 = var.sss(0), ss2 = var.sss(1);
  double v1 = var.svr(0), v2 = var.svr(1);
  double m1 = var.smean(0), m2 = var.smean(1);
  double md = m2 - m1;
  double vmd = (ss1 + ss2) / (double)df * (double)nt / (double)n1 / (double)n2;
  fout << "・順序効果：";
  put_mdci(fout, opt, md, df, vmd);
}

void crosovclass::put_dgd(ostream& fout, option& opt, int iidx, int ni, double mdi, double dij, double sed, double ci, double sedm, double dmin) const
// 処置差とその信頼区間と検出差を出力する
{
  fout << "  観測差=" << dij;
  if (opt.logt() == logt10) fout << "  対数変換時=log10(" << pow(10.0, dij);
  else fout << "  対数変換時=ln(" << exp(dij);
  fout << ")  " << opt.dg(iidx) << "に対する率=" << dij / mdi * 100.0 << "%\n";
  double dlow = dij - ci, dup = dij + ci;
  fout << "  " << opt.cc() * 100.0 << "%信頼区間:" << dij << "±" << ci << "(下限=" << dlow << "  上限=" << dup << ")\n";
  fout << "    対数変換時:下限=";
  if (opt.logt() == logt10) fout << "log10(" << pow(10.0, dlow) << ")  上限=log10(" << pow(10.0, dup) << ")\n";
  else fout << "ln(" << exp(dlow) << ")  上限=ln(" << exp(dup) << ")\n";
  fout << "    " << opt.dg(iidx) << "に対する率:下限=" << dlow / mdi * 100.0 << "%  上限=" << dup / mdi * 100.0 << "%\n";
  fout << "  " << ni << "例での検出差=" << dmin;
  if (opt.logt() == logt10) fout << "  対数変換時=log10(" << pow(10.0, fabs(dmin));
  else fout << "  対数変換時=ln(" << exp(dmin);
  fout << ")  " << opt.dg(iidx) << "に対する率=" << fabs(dmin / mdi * 100.0) << "%\n";
}

void crosovclass::put_logdn(ostream& fout, option& opt, double ta, double za, double zb, int ni, double dij, double sed, int dfr, double vr) const
// 対数検出差に関する検出力と必要例数を出力する
{
  double dd = opt.dl();
  double di = (za + zb) / dd;
  int n = (int)(di * di * vr + 1.9);
  double t1 = (dd / sed) - ta, t2 = 0.0;
  double pb = (t1 > 0.0) ? ttop(t1, dfr) * 0.5 : 1.0 - ttop(-t1, dfr) * 0.5;
  if (opt.logt() == logt10) fout << "  検出差log10(" << pow(10.0, dd);
  else fout << "  検出差ln(" << exp(dd);
  fout << ")=" << dd << " 1群=" << ni << "例の時の検出力=" << (1.0 - pb) * 100.0 << "%\n";
  if (opt.logt() == logt10) fout << "  検出差log10(" << pow(10.0, dd);
  else fout << "  検出差ln(" << exp(dd);
  fout << ")=" << dd << "の時の検定の1群必要例数=" << n << "\n";
  di = (ptoz(1.0 - opt.po()) + za) / dd;
  n = (int)(di * di * vr + 1.9);
  int df = 2 * n - 2;
  double ci = ptot(opt.al(), df) * sqrt(vr / (double)n);
  ldbl po0 = 0.0L;
  fout << "    観測差=0の時の絶対精度=" << ci << " 1群必要例数=" << n << "\n";
  double d = fabs(dij);
  if (d < dd)                  // 観測差が検出差よりも小さい時だけ出力
  {
    n = cal_nci(opt, ci, po0, dd, d, vr);
    fout << "    観測差=" << d << "の時の絶対精度2=" << ci << " 1群必要例数=" << n << "\n";
  }
}

void crosovclass::put_prodn(ostream& fout, option& opt, double ta, double za, double zb, int ni, double mdi, double dij, double sed, int dfr, double vr) const
// 割合検出差に関する検出力と必要例数を出力する
{
  double dd = fabs(mdi * opt.dp());
  double di = (za + zb) / dd;
  int n = int(di * di * vr + 1.9);
  double t1 = (dd / sed) - ta, t2 = 0.0;
  double pb = (t1 > 0.0) ? ttop(t1, dfr) * 0.5 : 1.0 - ttop(-t1, dfr) * 0.5;
  fout << "  検出差" << fabs(dd / mdi * 100.0) << "%=" << dd << " 1群=" << ni << "例の時の検出力=" << (1.0 - pb) * 100.0 << "%\n";
  fout << "  検出差" << fabs(dd / mdi * 100.0) << "%=" << dd << "の時の検定の1群必要例数=" << n << "\n";
  di = (ptoz(1.0 - opt.po()) + za) / dd;
  n = int(di * di * vr + 1.9);
  int df = 2 * n - 2;
  double ci = ptot(opt.al(), df) * sqrt(vr / double(n));
  ldbl po0 = 0.0L;
  fout << "    観測差=0の時の絶対精度=" << ci << " 1群必要例数=" << n << "\n";
  double d = fabs(dij);
  if (d < dd)                  // 観測差が検出差よりも小さい時だけ出力
  {
    n = cal_nci(opt, ci, po0, dd, d, vr);
    fout << "    観測差=" << d << "の時の絶対精度=" << ci << " 1群必要例数=" << n << "\n";
  }
}

int crosovclass::cal_nci(option& opt, double& ci, ldbl& po0, double dd, double d, double vr) const
// 信頼区間による必要例数を計算して返す
{
  if (dd < d) return 0;       // 計算不可能
  constexpr int Nmax = 10000; // 最大例数(繰り返し計算が終わらないため)
  double de = dd - d;
  double bci = 1.0 - opt.po() * de / dd;
  double di = (ptoz(opt.al()) + ptoz(bci)) / de;
  int n = int(di * di * vr + 0.9);
  int df = 2 * n - 2;
  double t0 = ptot(opt.al(), df);
  double sed = sqrt(vr / double(n));
  double tl = (-dd - d) / sed + t0, tu = (dd - d) / sed - t0;
  po0 = ncttolp(tu, df) - ncttolp(tl, df);
  while (po0 < (ldbl)(opt.po()) && n < Nmax)
  {
    n++;
    df = 2 * n - 2;
    t0 = ptot(opt.al(), df);  sed = sqrt(vr / double(n));
    tl = (-dd - d) / sed + t0;  tu = (dd - d) / sed - t0;
    po0 = ncttolp(tu, df) - ncttolp(tl, df);
  }
  ci = t0 * sed;                // 絶対精度
  return n;
}


void crosovclass::put_ttest(ostream& fout, option& opt, const varclass& var) const
// クロスオーバー用t検定結果を出力する
{
  int n1 = var.sn(0), n2 = var.sn(1);
  double ss1 = var.sss(0), ss2 = var.sss(1);
  double v1 = var.svr(0), v2 = var.svr(1);
  double m1 = var.smean(0), m2 = var.smean(1);
  fout << "\n○2標本t検定\n";
  string tstr("2時期の和の基礎統計量(順序効果)");
  put_tstat(fout, opt, tstr, n1, m1, ss1, v1, n2, m2, ss2, v2);
  n1 = var.dn(0);  n2 = var.dn(1);
  ss1 = var.dss(0);  ss2 = var.dss(1);
  v1 = var.dvr(0);  v2 = var.dvr(1);
  m1 = var.dmean(0);  m2 = -var.dmean(1);
  tstr = "2時期のクロス差の基礎統計量(時期効果)";
  fout << "\n";
  put_tstat(fout, opt, tstr, n1, m1, ss1, v1, n2, m2, ss2, v2);
  m1 = var.dmean(0);  m2 = var.dmean(1);
  tstr = "2時期の差の基礎統計量(処置効果)";
  fout << "\n";
  put_tstat(fout, opt, tstr, n1, m1, ss1, v1, n2, m2, ss2, v2);
}

void crosovclass::put_tstat(ostream& fout, option& opt, string& tstr, int n1, double m1, double ss1, double v1, int n2, double m2, double ss2, double v2) const
// 2標本t検定の基礎統計量と検定結果を出力する
{
  constexpr int LENline = 77;   // 区切り線の長さ
  constexpr int LENhed = 9;     // ヘッダの長さ
  fout << "群       :" << tstr << "\n";
  putline(fout, LENline);
  fout << left << setw(LENhed) << opt.g(0) << right << ":";
  putmsd(fout, n1, m1, v1);
  fout << left << setw(LENhed) << opt.g(1) << right << ":";
  putmsd(fout, n2, m2, v2);
  putline(fout, LENline);
  int nt = n1 + n2, df = n1 + n2 - 2;
  ldbl pr = 0.0L;
  bool welch = false;           // Welchの検定使用フラッグ
  string tname("等分散t検定");
  if (opt.tt() == testcor && n1 > 1 && n2 > 1) // Welchの検定使用
  { welch = true; }
  else if (opt.tt() == testaut) // 等分散性の両側F検定実行
  {
    double f = (v2 > 0.0) ? v1 / v2 : 0.0;
    int df1 = n1 - 1, df2 = n2 - 1;
    if (v1 < v2)
    {
      f = (v1 > 0.0) ? v2 / v1 : 0.0;
      df1 = n2 - 1;  df2 = n1 - 1;
    }
    pr = ftop(f, df1, df2);
    if (2.0L * pr <= (ldbl)(opt.al())) welch = true;
  }
  double md = m2 - m1;
  double vmd = (ss1 + ss2) / (double)df * (double)nt / (double)n1 / (double)n2;
  if (welch)
  {
    v1 /= (double)n1;  v2 /= (double)n2;
    vmd = v1 + v2;
    df = (int)(vmd * vmd / (v1 * v1 / (double)(n1 - 1) + v2 * v2 / (double)(n2 - 1)));
    tname = "不等分散t検定(Welchの検定)";
  }
  double t = md / sqrt(vmd);
  pr = ttop(t, df);
  fout << tname << ":t=" << t;
  fout << "  自由度=" << df << "  p=" << pr << pind(pr) << "\n";
  put_mdci(fout, opt, md, df, vmd);
}

void crosovclass::put_mdci(ostream& fout, option& opt, double md, int df, double vmd) const
// 平均値の差の信頼区間を出力する
{
  fout << "観測差の" << opt.cc() * 100.0 << "%信頼区間=";
  double ci = calmcci(mcnon, opt.cc(), 2, 1, df, 1, vmd) * 0.5;
  md *= 0.5;
  fout << md << "±" << ci << "(" << md - ci << "-" << md + ci << ")\n";
}

void crosovclass::put_utest(ostream& fout, option& opt, const varclass& var)
// クロスオーバー用U検定結果を出力する
{
  fout << "\n○ノンパラメトリック検定\n";
  ftable ftab(2, (var.sumset()).size()); // クロス前後の和と差の度数表
  sum_ftab(opt, ftab, sumorder, var);
  double u = 0.0, vu = 0.0, z = 0.0;
  testtype utype = utest(u, vu, z, ftab, opt.ut());
  ldbl pr = ztop(z);
  fout << "・順序効果:Wilcoxonの2標本検定(順位和検定、Mann-WhitneyのU検定";
  if (utype == testcor) fout << ",修正有";
  fout << ")\n";
  fout << "  正規分布z=" << z << "  有意確率p=" << pr << pind(pr) << "\n";
  put_umrdci(fout, opt, ftab, u, vu);
  ftab.resize(2, (var.difset()).size());
  sum_ftab(opt, ftab, sumperiod, var);
  utype = utest(u, vu, z, ftab, opt.ut());
  pr = ztop(z);
  fout << "\n・時期効果:Wilcoxonの2標本検定(順位和検定、Mann-WhitneyのU検定";
  if (utype == testcor) fout << ",修正有";
  fout << ")\n";
  fout << "  正規分布z=" << z << "  有意確率p=" << pr << pind(pr) << "\n";
  put_umrdci(fout, opt, ftab, u, vu);
  sum_ftab(opt, ftab, sumtreat, var);
  utype = utest(u, vu, z, ftab, opt.ut());
  pr = ztop(z);
  fout << "\n・処置効果:Wilcoxonの2標本検定(順位和検定、Mann-WhitneyのU検定";
  if (utype == testcor) fout << ",修正有";
  fout << ")\n";
  fout << "  正規分布z=" << z << "  有意確率p=" << pr << pind(pr) << "\n";
  put_umrdci(fout, opt, ftab, u, vu);
}

void crosovclass::put_umrdci(ostream& fout, option& opt, ftable& ftab, double u, double vu) const
// 順位平均値の差の信頼区間とU値の信頼区間を出力する
{
  int n1 = ftab.frtotal(1), n2 = ftab.frtotal(2);
  int nt = n1 + n2;
  double ta = u + (double)n2 * (double)(n2 + 1) * 0.5;
  double r2 = ta / (double)n2, tt = (double)nt * (double)(nt + 1) * 0.5;
  double r1 = (tt - ta) / (double)n1;
  double drm = r2 - r1;
  if (fabs(drm) < FDeps) drm = 0.0;
  double vd = vu / (double)n1 / (double)n2 * (double)nt * (1.0 / (double)n1 + 1.0 / (double)n2);
  double ci = calmcci(mcnon, opt.cc(), ftab.nr(), 1, 0, 1, vd);
  double dl = drm - ci, du = drm + ci;
  double pro = 100.0 / (double)nt;
  fout << "  順位平均値の差の" << opt.cc() * 100.0 << "%信頼区間(confidence interval)\n";
  fout << "    " << drm << "(" << drm * pro << "%)±" << ci;
  fout << "(" << dl << "(" << dl * pro << "%)-" << du << "(" << du * pro << "%))\n";
  ci = calmcci(mcfis, opt.cc(), ftab.nr(), 1, 0, 1, vu);
  dl = u - ci;  du = u + ci;
  if (dl < 0.0) dl = 0.0;
  if (du > (double)n1 * (double)n2) du = (double)n1 * (double)n2;
  pro = 100.0 / (double)n1 / (double)n2;
  fout << "  U値の" << opt.cc() * 100.0 << "%信頼区間(confidence interval)\n";
  fout << "    " << u << "(" << u * pro << "%)±" << ci;
  fout << "(" << dl << "(" << dl * pro << "%)-" << du << "(" << du * pro << "%))\n";
}

void crosovclass::sum_ftab(option& opt, ftable& ftab, sumtype stype, const varclass& var)
// クロスオーバー用U検定用用にクロス前後の和／差／クロス差を度数集計する
{
  ftab.fclear();
  for (auto& cdatai : datalst)
  {
    int rno = opt.search_g(cdatai.sdata(0)) + 1;
    if (rno == 1 || rno == 2)   // 群1または2だけ集計
    {
      double fdpre = cdatai.fdata(1), fdpost = cdatai.fdata(2);
      if (fdpre != BLANKf && fdpost != BLANKf)
      {
        double fdata = 0.0;
        int cno = 0;
        switch (stype)
        {
        case sumorder : fdata = fdpre + fdpost;
          cno = search_data(var.sumset(), fdata);
          break;
        case sumperiod : fdata = (rno == 1) ? fdpost - fdpre : fdpre - fdpost;
          cno = search_data(var.difset(), fdata);
          break;
        case sumtreat : fdata = fdpost - fdpre;
          cno = search_data(var.difset(), fdata);
        }
        if (cno > 0) ftab.f(rno, cno, ftab.f(rno, cno) + 1);
      }
    }
  }
  ftab.set();
}

int crosovclass::search_data(const set<double>& targset, double fdata) const
// 数字データをサーチして番号を返す
{
  int sno = 1;
  bool hit = false;
  set<double>::const_iterator sitr = targset.begin();
  while (hit == false && sitr != targset.end())
  {
    if (sgndiff(fdata, *sitr) == 0) hit = true; // 誤差範囲なら一致
    else { ++sno;  ++sitr; }
  }
  return hit ? sno : -1;        // -1:なし
}

int crosovclass::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() - 1) != 0) throw option_error("実行時オプションエラー！");
  return nc;
}

void crosovclass::put(option& opt) // 集計結果を出力する
{
  if (opt.res("res") == false) return;
  put_resid(cout);
  varclass var(id.ni() - 1);        // 計算用変数クラス
  var.cal(opt, datalst);
  if (var.tn() < 1)
  {
    cout << "※有効なデータなし！\n\n" << endl;
    return;
  }
  if (opt.res("stat")) put_stat(cout, opt, var);
  if (var.resvar() == reserr)
  { cout << "※分散分析不可能！\n"; }
  else
  {
    if (opt.res("at")) put_coatab(cout, var);
    put_dgest(cout, opt, var);
    if (opt.res("tt")) put_ttest(cout, opt, var);
    if (opt.res("ut")) put_utest(cout, opt, var);
  }
  cout << "\n" << endl;
}

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