// fissim.C  Version 7.42  Copyleft (c) 杉本典夫
// Time-stamp: <2025-11-03 11:24:31 tac>
//
// Fisherの正確検定のシミュレーション
//
// 使用ライブラリ:libdans.so, libfun.so
// コンパイル法  :g++ -O -s -W fissim fissim.C -I/usr/local/include
//               -L/usr/local/lib -ldans -lfun

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

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

const vecstr USAGE =            // プログラム使用法
{
  "fissim  Version 7.42         Copyleft (Vo\\oV) 2025 杉本典夫",
  "書式:fissim  [<オプション> …]  [<データファイル> …]",
  "機能:Fisherの正確検定のシミュレーション。",
  "<オプション>",
  "  -cc[=<信頼係数>]:信頼係数を指定する。",
  "    (0≦信頼係数＜1,オプション省略時:-cc=0.95)",
  "  -f=<ファイル名>:指定されたファイルからオプションを読み込む。",
  "  -h[=<コマンド名>]:オンラインマニュアルを表示する。(-h:-h=man)",
  "  -n[{1|2}][=<例数>]:群1と群2の例数を指定する。",
  "    (0＜例数,オプション省略時:-n=10)",
  "  -nc1[=<例数>]:分類1の合計例数を指定する。",
  "    (0＜例数,オプション省略時:合計例数の半分)",
  "  -two[={-|+}]:両側検定の指定。(-two:two=+)",
  "    -:片側検定  +:両側検定(オプション省略時)"
};

const string Title = "=== Fisherの正確検定のシミュレーション ==="; // 表題
constexpr int Ndef = 10;        // 省略時例数

double cal_dpci(double& dpci, double p1, double p2, int n1, int n2, double za, bool cor)
// 正規近似法による信頼区間幅と出現率差を返す
{
  if (n1 < 1 || n2 < 1 || za <= 0.0) return 0.0;
  int nt = n1 + n2;
  double v = p1 * (1.0 - p1) / (double)n1 + p2 * (1.0 - p2) / (double)n2;
  dpci = za * sqrt(v);
  if (cor && sgn(dpci) != 0) dpci += 0.5 / (double)n1 + 0.5 / (double)n2;
  return p1 - p2;               // 出現率差を返す
}

class option                    // オプションクラス
{
private:
  string v_path;                // プログラムのパス
  string v_pname;               // プログラム名
  double v_cc;                  // 信頼係数
  int v_n1;                     // 群1の例数
  int v_n2;                     // 群2の例数
  int v_nc1;                    // 分類1の合計例数
  bool v_two;                   // 両側検定フラッグ
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; }
  int n1() const                // 群1の例数を返す
  { return v_n1; }
  int n2() const                // 群2の例数を返す
  { return v_n2; }
  int nc1() const               // 分類1の合計例数を返す
  { return v_nc1; }
  bool two() const              // 両側検定フラッグを返す
  { return v_two; }
  void setdef()                 // 省略時オプション指定を設定する
  { getoptdef(*this); }
  void set(const string& optstr, char devc = OPTdev);
  // オプション指定によってオプションを設定する
  int check();                  // オプションのチェック
};

option::option(const char* c_pname) : v_path(), v_pname(), v_cc(CCdef), v_n1(Ndef), v_n2(Ndef), v_nc1(0), v_two(true)
// コンストラクタ
{
  if (c_pname == nullptr) return; // 初期設定せず
  v_pname = getname(c_pname);     // パス名なしのプログラム名取得
  v_path = c_pname;
  strsize_t idx = v_path.find(v_pname);
  if (idx != string::npos) v_path.erase(idx); // プログラム名削除
}

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] == "n")       // 群1=群2の例数の指定
  {
    char* ptr = nullptr;
    v_n1 = opt[1].empty() ? Ndef : strtol(opt[1].c_str(), &ptr, 0);
    if (v_n1 < 1) err = badopt(pname(), optstr, USAGE, 1.0, (double)INT_MAX);
    v_n2 = v_n1;
  }
  else if (opt[0] == "n1")      // 群1の例数の指定
  {
    char* ptr = nullptr;
    v_n1 = opt[1].empty() ? Ndef : strtol(opt[1].c_str(), &ptr, 0);
    if (v_n1 < 1) err = badopt(pname(), optstr, USAGE, 1.0, (double)INT_MAX);
  }
  else if (opt[0] == "n2")      // 群2の例数の指定
  {
    char* ptr = 0;
    v_n2 = opt[1].empty() ? Ndef : strtol(opt[1].c_str(), &ptr, 0);
    if (v_n2 < 1) err = badopt(pname(), optstr, USAGE, 1.0, (double)INT_MAX);
  }
  else if (opt[0] == "nc1")     // 分類1の合計例数の指定
  {
    char* ptr = nullptr;
    v_nc1 = opt[1].empty() ? Ndef : strtol(opt[1].c_str(), &ptr, 0);
    if (v_nc1 < 1) err = badopt(pname(), optstr, USAGE, 1.0, (double)INT_MAX);
  }
  else if (opt[0] == "two")     // 両側検定の指定
  {
    if (opt[1].empty() || opt[1] == "+") v_two = true;
    else if (opt[1] == "-") v_two = false;
    else err = badopt(pname(), optstr, USAGE);
  }
  else                          // 間違ったオプション指定
  { err = nonopt(pname(), optstr, USAGE); }
  if (err > 0) throw help_except("ヘルプ表示");
  if (err < 0) throw option_error(string("オプション指定違反→") + optstr);
}

int option::check()             // オプションのチェック
{
  if (v_nc1 >= v_n1 + v_n2) return disperr(pname(), "分類1の合計例数≧合計例数違反！");
  return 0;                     // エラーなし
}

class varclass                  // 変数クラス
{
private:
  int v_n1;                     // 群1の例数
  int v_n2;                     // 群2の例数
  int v_nc1;                    // 分類1の例数
  int v_nc2;                    // 分類2の例数
  int v_nmin;                   // 周辺度数最小値
  int v_nlow;                   // 下限棄却域度数
  int v_nup;                    // 上限棄却域度数
  double v_dp;                  // 母出現率差
  double v_dpl;                 // 母出現率差下限
  double v_dpu;                 // 母出現率差上限
  double v_pl;                  // 上側棄却域の累積確率
  double v_pu;                  // 下側棄却域の累積確率
  vecdbl fisvec;                // 超幾何分布ベクトル
  vecdbl v_ndpvec; // 正規近似法(全体の出現率使用・連続修正有)のクロス集計表
  vecdbl v_ndcvec;             // 正規近似法(連続修正有)のクロス集計表
  vecdbl v_ndvec;              // 正規近似法(連続修正無)のクロス集計表
public:
  varclass(int n1, int n2, int nc1, int nc2, double cc0, bool two);
  // コンストラクタ
  ~varclass()                   // デストラクタ
  {}
  void set(double cc0, bool two); // 変数をセットする
  int n1() const                  // 群1の例数を返す
  { return v_n1; }
  int n2() const                // 群2の例数を返す
  { return v_n2; }
  int nc1() const               // 分類1の例数を返す
  { return v_nc1; }
  int nc2() const               // 分類2の例数を返す
  { return v_nc2; }
  int nmin() const              // 周辺度数最小値を返す
  { return v_nmin; }
  int nlow() const              // 下限棄却域度数を返す
  { return v_nlow; }
  int nup() const               // 上限棄却域度数を返す
  { return v_nup; }
  double dp() const             // 母出現率差を返す
  { return v_dp; }
  double dpl() const            // 母出現率差下限を返す
  { return v_dpl; }
  double dpu() const            // 母出現率差上限を返す
  { return v_dpu; }
  double pl() const             // 下側棄却域の累積確率を返す
  { return v_pl; }
  double pu() const             // 上側棄却域の累積確率を返す
  { return v_pu; }
  const vecdbl& ndpvec() const
  // 正規近似法(全体の出現率使用・連続修正有)のクロス集計表を返す
  { return v_ndpvec; }
  const vecdbl& ndcvec() const // 正規近似法(連続修正有)のクロス集計表を返す
  { return v_ndcvec; }
  const vecdbl& ndvec() const // 正規近似法(連続修正無)のクロス集計表を返す
  { return v_ndvec; }
};

varclass::varclass(int n1, int n2, int nc1, int nc2, double cc0, bool two) : v_n1(0), v_n2(0), v_nmin(0), v_nc1(0), v_nc2(0), v_nlow(0), v_nup(0), v_dp(0.0), v_dpl(0.0), v_dpu(0.0), v_pl(0.0), v_pu(0.0), fisvec(), v_ndpvec(4, 0.0), v_ndcvec(4, 0.0), v_ndvec(4, 0.0)
// コンストラクタ
{
  int nt = n1 + n2;
  v_n1 = min(n1, n2);           // vn1:例数が少ない方
  v_n2 = nt - v_n1;
  v_nc1 = min(nc1, nc2);        // vnc1:例数が少ない方
  v_nc2 = nt - v_nc1;
  v_nmin = min(v_n1, v_nc1);    // 周辺度数最小値
  fisvec.resize(v_nmin + 1, 0.0);
  double alow = two ? (1.0 - cc0) * 0.5 : 0.0;
  double aup = cc0 + alow;
  double cipr = 0.0;
  bool lowsch = two ? true : false;
  bool upsch = true;
  v_nlow = 0;  v_nup = v_nmin;  v_pl = 0.0;
  for (int i = 0; i <= v_nmin; ++i)
  {
    int n11 = i, n12 = v_n1 - i;
    int n21 = v_nc1 - n11, n22 = v_nc2 - n12;
    double fpr = fisden(n11, n12, n21, n22);
    fisvec[i] = fpr;
    if (lowsch)                 // 信頼限界下限探索
    {
      if (v_pl + fpr <= alow) v_pl += fpr;
      else { v_nlow = i;  lowsch = false; }
    }
    if (upsch)                  // 信頼限界上限探索
    {
      cipr += fpr;
      if (cipr >= aup) { v_nup = i;  upsch = false; }
    }
  }
  double dp1 = (double)v_nlow / (double)v_n1;
  int n21 = v_nc1 - v_nlow;
  double dp2 = (double)n21 / (double)v_n2;
  v_dpl = dp1 - dp2;            // 出現率差下限
  dp1 = (double)v_nup / (double)v_n1;
  n21 = v_nc1 - v_nup;
  dp2 = (double)n21 / (double)v_n2;
  v_dpu = dp1 - dp2;                       // 出現率差上限
  v_pu = (v_nup < v_nmin) ? 1.0 - cipr : 0.0; // 上側棄却域の累積確率
  v_dp = 0.0;                                 // 母出現率差
}

void varclass::set(double cc0, bool two) // 変数をセットする
{
  int nt = n1() + n2();
  double p0 = (double)nc1() / (double)nt;
  double alpha = 1.0 - cc0;
  if (two == false) alpha *= 2.0;
  double za = ptoz(alpha);
  size_t nocip = 0, nocis = 0;
  for (int i = 0; i <= nmin(); ++i)
  {
    int n11 = i, n12 = n1() - i;
    int n21 = nc1() - n11, n22 = nc2() - n12;
    double p1 = (double)n11 / (double)n1();
    double p2 = (double)n21 / (double)n2();
    double dpi = p1 - p2, dpci = 0.0;
    if (two) nocip = (dpl() <= dpi && dpi <= dpu()) ? 0 : 1;
    else nocip = (dpi <= dpu()) ? 0 : 1; // 母集団の信頼区間チェック
    cal_dpci(dpci, p0, p0, n1(), n2(), za, true); // 母出現率を使用した正規近似法
    double dpli = dpi - dpci, dpui = dpi + dpci;
    if (two) nocis = (dpli <= dp() && dp() <= dpui) ? 0 : 1;
    else nocis = (dpli <= dp()) ? 0 : 1;
    size_t idx = nocip * 2 + nocis;
    v_ndpvec[idx] += fisvec[i];
    dpi = cal_dpci(dpci, p1, p2, n1(), n2(), za, true); // 正規近似法(連続修正有)
    dpli = dpi - dpci;  dpui = dpi + dpci;
    if (two) nocis = (dpli <= dp() && dp() <= dpui) ? 0 : 1;
    else nocis = (dpli <= dp()) ? 0 : 1;
    idx = nocip * 2 + nocis;
    v_ndcvec[idx] += fisvec[i];
    dpi = cal_dpci(dpci, p1, p2, n1(), n2(), za, false); // 正規近似法(連続修正無)
    dpli = dpi - dpci;  dpui = dpi + dpci;
    if (two) nocis = (dpli <= dp() && dp() <= dpui) ? 0 : 1;
    else nocis = (dpli <= dp()) ? 0 : 1;
    idx = nocip * 2 + nocis;
    v_ndvec[idx] += fisvec[i];
  }
}

class fissimclass               // fissimメインクラス
{
private:
  void put_fis(ostream& fout, const option& opt, const varclass& var) const;
  // 検定結果と推定結果のクロス集計表を出力する
  void put_table(ostream& fout, const option& opt, const vecdbl& tmpvec) const;
  // 分割表と診断率を出力する
public:
  fissimclass()                 // コンストラクタ
  {}
  ~fissimclass()                // デストラクタ
  {}
  void put(option& opt);        // シミュレーション結果を出力する
};

void fissimclass::put_fis(ostream& fout, const option& opt, const varclass& var) const
// 検定結果と推定結果のクロス集計表を出力する
{
  puttitle(fout, Title);
  fout << "\n○検定:Fisherの正確検定  群1の例数=" << opt.n1() << "  群2の例数=" << opt.n2() << "\n";
  double pro = opt.cc() * 100.0;
  int nt = opt.n1() + opt.n2();
  int nc1 = (opt.nc1() > 0) ? opt.nc1() : nt / 2; // nc1はオプションで指定した値
  int nc2 = nt - nc1;
  fout << "  分類1の合計例数=" << nc1 << "  分類2の合計例数=" << nc2 << "  母出現率差=" << var.dp();
  if (opt.two()) fout << "  両側検定\n";
  else fout << "  片側検定\n";
  int n21 = var.nc1() - var.nlow();
  fout << "  " << pro << "%信頼区間=" << var.dpl() << "(" << var.nlow() << "/" << var.n1() << "-" << n21 << "/" << var.n2() << ")";
  n21 = var.nc1() - var.nup();
  fout << "〜" << var.dpu() << "(" << var.nup() << "/" << var.n1() << "-" << n21 << "/" << var.n2() << ")\n";
  fout << "  " << 100.0 - pro << "%棄却域:下側累積確率=" << var.pl() << "  上側累積確率=" << var.pu() << "\n";
  double p0 = (double)var.nc1() / (double)nt;
  double alpha = 1.0 - opt.cc();
  if (opt.two() == false) alpha *= 2.0;
  double cc = 1.0 - alpha, za = ptoz(alpha);
  double dpci = 0.0;
  cal_dpci(dpci, p0, p0, var.n1(), var.n2(), za, true);
  double dpli = var.dp() - dpci, dpui = var.dp() + dpci;
  if (opt.two() == false) dpui = 1.0;
  fout << "\n・推定:正規近似法(母出現率使用)\n";
  fout << "  出現率差=" << var.dp() << "の時の" << pro << "%信頼区間=" << dpli << "〜" << dpui << "\n";
  put_table(fout, opt, var.ndpvec());
  int n11 = (double)var.n1() * p0;
  n21 = var.nc1() - n11;
  double p1 = (double)n11 / (double)var.n1();
  double p2 = (double)n21 / (double)var.n2();
  double dpi = cal_dpci(dpci, p1, p2, var.n1(), var.n2(), za, true);
  dpli = dpi - dpci;  dpui = dpi + dpci;
  if (opt.two() == false) dpui = 1.0;
  fout << "\n・推定:正規近似法(連続修正有)\n";
  fout << "  出現率差=" << dpi << "の時の" << pro << "%信頼区間=" << dpli << "〜" << dpui << "\n";
  put_table(fout, opt, var.ndcvec());
  dpi = cal_dpci(dpci, p1, p2, var.n1(), var.n2(), za, false);
  dpli = dpi - dpci;  dpui = dpi + dpci;
  if (opt.two() == false) dpui = 1.0;
  fout << "\n・推定:正規近似法(連続修正無)\n";
  fout << "  出現率差=" << dpi << "の時の" << pro << "%信頼区間=" << dpli << "〜" << dpui << "\n";
  put_table(fout, opt, var.ndvec());
}

void fissimclass::put_table(ostream& fout, const option& opt, const vecdbl& tmpvec) const
// 分割表と診断率を出力する
{
  double r1total = tmpvec[0] + tmpvec[1], r2total = tmpvec[2] + tmpvec[3];
  double c1total = tmpvec[0] + tmpvec[2], c2total = tmpvec[1] + tmpvec[3];
  double pro = opt.cc() * 100.0;
  fout << "検定結果＼推定結果    信頼区間      棄却域        合計\n";
  int wlen = 54, w = 12;
  putline(fout, wlen);
  fout << "信頼区間          " << setw(w) << tmpvec[0] << setw(w) << tmpvec[1] << setw(w) << r1total << "\n";
  fout << "棄却域            " << setw(w) << tmpvec[2] << setw(w) << tmpvec[3] << setw(w) << r2total << "\n";
  putline(fout, wlen);
  fout << "合計              " << setw(w) << c1total << setw(w) << c2total << setw(w) << 1 << "\n";
  pro = (r1total > 0.0) ? tmpvec[0] / r1total : 0.0;
  fout << "感度=" << tmpvec[0] << "/" << r1total << "=" << pro;
  pro = (r2total > 0.0) ? tmpvec[3] / r2total : 0.0;
  fout << " 特異度=" << tmpvec[3] << "/" << r2total << "=" << pro << "\n";
  pro = tmpvec[0] + tmpvec[3];
  fout << "正診率=" << tmpvec[0] << "+" << tmpvec[3] << "=" << pro << "\n";
}

void fissimclass::put(option& opt) // シミュレーション結果を出力する
{
  if (opt.check() != 0) throw option_error("実行時オプションエラー！");
  int nt = opt.n1() + opt.n2();
  int nc1 = (opt.nc1() > 0) ? opt.nc1() : nt / 2;
  int nc2 = nt - nc1;
  varclass var(opt.n1(), opt.n2(), nc1, nc2, opt.cc(), opt.two());
  // 変数オブジェクト
  var.set(opt.cc(), opt.two()); // 変数をセットする
  put_fis(cout, opt, var);
  cout << "\n" << endl;
}

int main(int argc, char* argv[]) // メインルーチン
{
  fissimclass fissim;           // メインクラス変数
  option opt(argv[0]);          // オプション
  try                           // 例外監視
  {
    opt.setdef();               // オプションの初期設定
    bool fok = true;            // ファイル処理モードフラッグ
    for (int i = 1; i < argc; ++i)
    {
      if (*(argv[i]) == CHRopt) opt.set(argv[i] + 1);
    }
    fissim.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;
}
