// dpplot.C Version 7.44  Copyleft (c) 杉本典夫
// Time-stamp: <2025-10-29 09:31:29 tac>
//
// 群項目と集計項目からDP-plot分析を行う。
//
// 使用ライブラリ:libdans.so, libfun.so
// コンパイル法 :g++ -O -s -W dpplot dpplot.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 =            // プログラム使用法表示用
{
  "dpplot  Version 7.44         Copyleft(Vo\\oV) 2025 杉本典夫",
  "書式:dpplot  [<オプション> …]  [<データファイル> …]",
  "機能:最初の項目を群項目、次の項目を集計項目としてDP-plot分析を行う",
  "<オプション>",
  "  -cc[=<信頼係数>]:信頼係数を指定する。(0≦信頼係数<1)",
  "    -cc:-cc=0.95(オプション省略時)",
  "  -co=<境界値>:指定した境界値の診断率を出力する。",
  "  -f=<ファイル名>:指定されたファイルからオプションを読み込む。",
  "  -g{1|2}[=<群コード>]:群1(疾患群)または群2(正常群)のコードを指定する。",
  "  -h[=<コマンド名>]:オンラインマニュアルを表示する。(-h:-h=man)",
  "  -lim={auto|low|up}:境界値の上限／下限を指定する。",
  "    auto:自動判定(オプション省略時)  low:下限  up:上限",
  "  -ngp={ac|dpp|log|pnf|pv|snp}[,…]:Ngraph用出力データの指定。",
  "    ac:正診率曲線  dpp:常用DP-plot用各種曲線  log:ロジスティック曲線",
  "    pnf:陽性率/陰性率曲線  pv:PPV/NPV曲線  snp:感度/特異度曲線",
  "    オプション省略時:全て指定",
  "  -ngpname[=<出力ファイル名>]:Ngraph用データを出力するファイル名を指定する。",
  "    オプション省略時:データ出力せず  -ngpname:-ngpname=dpplot_ngp.dat",
  "  -pi[=<群1の事前確率>]:群1の事前確率を指定します。(0<事前確率<1)",
  "  -res=[-][{def|all|ac[<正診率>]|dpa|stat|uda}][,…]",
  "      :解析結果の出力項目を指定する。",
  "    -:-を付けた項目を出力しない(-res=-:全ての項目を出力しない)",
  "    def:オプション省略時と同じ項目を出力する  all:全ての項目を出力する",
  "    ac[<正診率>]:指定した正診率になる時の境界値",
  "    dpa:DP-plot分析  stat:基礎統計量  uda:単変量線形判別分析",
  "    オプション省略時:データに合わせて自動判定",
  "  -xrange=<最小値>-<最小値>:X軸の最小値と最大値を指定する。"
};

const string Title = "=== DP-plot(Diagnostic Performance plot)分析 ==="; // 表題
const string NGPNAMEdef = "dpplot_ngp.dat"; // 省略時Ngraphファイル名
constexpr uint RDdef = 50;      // 座標軸の分割数
constexpr size_t Nco = 9;       // 境界値の個数
constexpr double Rerror = 1.0e-8; // 収束限界
constexpr double ACdef = 0.7;     // 省略時正診率
constexpr double PIdef = 0.5;     // 省略時事前確率

enum limtype                    // 境界値の種類
{
  limnon,                       // 指定なし
  limauto,                      // 自動判定
  limlow,                       // 下限
  limup,                        // 上限
  limetc                        // その他
};

class option                    // オプション集合体
{
private:
  string v_path;                // プログラムのパス
  string v_pname;               // プログラム名
  double v_ac;                  // 正診率
  double v_cc;                  // 信頼係数
  double v_co;                  // 境界値
  vecstr gvec;                  // 群コードベクトル
  limtype v_lim;                // 境界値の種類
  string v_ngpname;             // Ngraphファイル名
  lstr ngplst;                  // Ngraph出力リスト
  mapsb ngpmap;                 // Ngraph出力の有無マップ
  double v_pi;                  // 群1の事前確率
  lstr reslst;                  // 分析結果出力リスト
  mapsb resmap;                 // 分析結果出力の有無マップ
  pairdbl xrange;               // X軸の最小値と最大値
public:
  option(const char* c_pname = nullptr); // コンストラクタ
  ~option()                              // デストラクタ
  {}
  const string& path() const    // プログラムのパスを返す
  { return v_path; }
  const string& pname() const   // プログラム名を返す
  { return v_pname; }
  double ac() const             // 正診率を返す
  { return v_ac; }
  double cc() const             // 信頼係数を返す
  { return v_cc; }
  double co() const             // 境界値を返す
  { return v_co; }
  const string& g(int idx) const // 群コードを返す
  { return gvec[idx]; }
  limtype lim() const           // 境界値の種類を返す
  { return v_lim; }
  const string& ngpname() const // Ngraphデータファイル名を返す
  { return v_ngpname; }
  bool ngp(const string& code); // Ngraph出力フラッグを返す
  double pi() const             // 群1の事前確率を返す
  { return v_pi; }
  bool res(const string& code); // 分析結果出力フラッグを返す
  double xmin() const           // X軸の最小値を返す
  { return xrange.first; }
  double xmax() const           // X軸の最大値を返す
  { return xrange.second; }
  void setdef()                 // 省略時オプション指定を設定する
  { getoptdef(*this); }
  void set(const string& optstr, char devc = OPTdev);
  // オプション指定によってオプションを設定する
  int search_g(const string& gcode) const; // 群コードを検索する
  int check(size_t ni);         // オプションのチェックとセット
  int set_ngp();                // -ngpオプションのセット
  int set_res();                // -resオプションのセット
};

option::option(const char* c_pname) : v_path(), v_pname(), v_ac(ACdef), v_cc(CCdef), v_co(BLANKf), gvec(2), v_lim(limauto), v_ngpname(), ngplst(), ngpmap(), v_pi(BLANKf), reslst(), resmap(), xrange(BLANKf, BLANKf)
// コンストラクタ
{
  if (c_pname == 0) return;     // 初期設定せず
  gvec[0] = "1";               // 群1のコードだけ初期設定
  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::ngp(const string& code) // Ngraph出力フラッグを返す
{
  mapsb_citr ngppos = ngpmap.find(code);
  if (ngppos == ngpmap.end()) return false; // code無し！
  return ngpmap[code];                     // true/falseを返す
}

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] == "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] == "co")      // 境界値の指定
  {
    v_co = opt[1].empty() ? BLANKf : strtodf(opt[1]);
    if (v_co == BLANKf) err = badopt(pname(), optstr, USAGE, -FDmax, FDmax);
  }
  else if (opt[0] == "f")       // オプション指定ファイル
  { err = fgetopt(*this, opt[1], true); }
  else if (opt[0] == "g1")      // 群1のコード指定
  {
    if (opt[1].empty()) gvec[0].erase();
    else gvec[0] = opt[1];
  }
  else if (opt[0] == "g2")      // 群2のコード指定
  {
    if (opt[1].empty()) gvec[1].erase();
    else gvec[1] = opt[1];
  }
  else if (opt[0] == "h")       // オンラインヘルプ表示
  { err = help(pname(), USAGE, opt[1]); }
  else if (opt[0] == "lim")
  {
    if (opt[1].empty() || opt[1] == "auto") v_lim = limauto;
    else if (opt[1] == "low") v_lim = limlow;
    else if (opt[1] == "up") v_lim = limup;
    else err = badopt(pname(), optstr, USAGE);
  }
  else if (opt[0] == "ngp")     // Ngraph出力項目指定
  {
    if (strsplit(ngplst, opt[1], CHRdelm) < 1) err = badopt(pname(), optstr, USAGE);
  }
  else if (opt[0] == "ngpname")     // Ngraphファイル名指定
  {
    if (opt[1].empty()) v_ngpname = NGPNAMEdef;
    else v_ngpname = opt[1];
    if (v_ngpname.size() < 3) err = badopt(pname(), optstr, USAGE);
  }
  else if (opt[0] == "pi")      // 群1の事前確率指定
  {
    v_pi = opt[1].empty() ? PIdef : strtodf(opt[1]);
    if (v_pi <= 0.0 || 1.0 <= v_pi) 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 if (opt[0] == "xrange")  // X軸の最小値と最大値の指定
  {
    if (getrange(xrange.first, xrange.second, opt[1]) == 0) // 範囲指定無
    { err = badopt(pname(), optstr, USAGE, -FDmax, FDmax); }
    else if (xrange.first >= xrange.second) // 範囲指定不適切
    { err = badopt(pname(), optstr, USAGE, -FDmax, FDmax); }
  }
  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_ngp() < 0) return -2; // -ngpオプションのセット
  if (set_res() < 0) return -3; // -resオプションのセット
  return 0;                     // エラーなし
}

int option::set_ngp()           // -ngpオプションのセット
{
  ngpmap["ac"] = ngpmap["dpp"] = ngpmap["log"] = ngpmap["pnf"] = ngpmap["pv"] = ngpmap["snp"] = true;
  if (ngplst.size() > 0)
  {
    if (v_ngpname == "") return disperr(pname(), "-ngpnameオプションの指定が必要！");
    ngpmap["ac"] = ngpmap["dpp"] = ngpmap["log"] = ngpmap["pnf"] = ngpmap["pv"] = ngpmap["snp"] = false;
    for (const string& ngpstr : ngplst)
    {
      mapsb_citr ngppos = ngpmap.find(ngpstr);
      if (ngppos == ngpmap.end())
      { return disperr(pname(), "-ngp=" + ngpstr + "オプション指定が不適当！"); }
      else
      { ngpmap[ngpstr] = true; }
    }
  }
  return 0;                     // エラーなし
}

int option::set_res()           // -resオプションのセット
{
  resmap["res"] = resmap["def"] = resmap["dpa"] = resmap["stat"] = resmap["uda"] = true;
  resmap["ac"] = false;
  if (reslst.size() > 0)
  {
    resmap["def"] = resmap["dpa"] = resmap["stat"] = resmap["uda"] = 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["dpa"] = resmap["stat"] = resmap["uda"] = true;
      }
      else if (resstr == "all")
      {
        resmap["def"] = resmap["ac"] = resmap["dpa"] = resmap["stat"] = resmap["uda"] = true;
      }
      else if (resstr[0] == 'a' && resstr[1] == 'c' && resstr.size() > 2)
      {
        v_ac = strtodf(resstr.substr(2)); // 正診率
        if (v_ac < 0.0 || 1.0 < v_ac) return disperr(pname(), "-res=" + resstr + "オプション指定が不適当！");
        resmap["ac"] = true;
      }
      else if (respos == resmap.end())
      { return disperr(pname(), "-res=" + resstr + "オプション指定が不適当！"); }
      else
      { resmap[resstr] = true; }
    }
  }
  return 0;                     // エラーなし
}

int option::search_g(const string& gcode) const // 群コードを検索する
{
  if (gcode == BLANKs) return -1; // 欠損データはヒットさせない
  int idx = -1;
  if (gvec[0].empty())         // 群1のコード指定がない場合
  {
    idx = 0;
    if (gcode == gvec[1]) idx = 1;
  }
  else if (gcode == gvec[0])
  { idx = 0; }
  else if (gvec[1].empty())    // 群2のコード指定がない場合
  { idx = 1; }
  else if (gcode == gvec[1])
  { idx = 1; }
  return idx;                   // -1:ヒットせず,0:群1,1:群2
}

class data_less                 // ケースデータの比較用関数オブジェクト
{
private:
  int kidx;                     // ソートキー項目のインデックス
public:
  data_less(int kidx0) : kidx(kidx0) // コンストラクタ
  {}
  ~data_less()                  // デストラクタ
  {}
  bool operator()(datac& cdata1, datac& cdata2) const; // ケースデータを比較する
};

bool data_less::operator()(datac& cdata1, datac& cdata2) const
// ケースデータを比較する
{
  int comp = 0;
  double fdata1 = cdata1.fdata(kidx);
  double fdata2 = cdata2.fdata(kidx);
  if (fdata1 < fdata2) comp = -1;
  else if (fdata1 > fdata2) comp = 1;
  return (comp < 0) ? true : false; // 小さいときだけtrueを返す
}

class diag                      // 診断率クラス
{
private:
  double v_pi;                   // 群1の事前確率
  bool v_uporder;                // データの昇順/下順の区別
  vecdbl covec;                  // 境界値
  vecdbl snvec;                  // 感度
  vecdbl spvec;                  // 特異度
  vecdbl acvec;                  // 正診率
  vecdbl ppvvec;                 // PPV
  vecdbl npvvec;                 // NPV
  double v_d2;                   // マハラノビスの汎距離
  double v_ppr;                  // 誤判別確率
  double v_aucn;                 // 理論的AUC
  restype v_resdiag;             // 計算結果の種類
  vector<statv<double>> statvec; // 基礎統計量ベクトル
  void search_acc(option& opt, ldata& datalst); // 実測値正診率点の探索
  void cal_uda(option& opt);                    // 単変量線形判別分析
  double cal_cospmin(double m1, double sd1, double m2, double sd2);
  // 理論的|SN-SP|≒0の時の境界値を計算する
  double cal_coacmax(double m1, double sd1, double m2, double sd2);
  // 理論的正診率最大値の時の境界値を計算する
  double cal_copnmin(double m1, double sd1, double m2, double sd2);
  // 理論的|PPV-NPV|≒0の時の境界値を計算する
  int search_coacc(vecdbl& covec, double acc, double m1, double sd1, double m2, double sd2);
  // 理論的正診率点の探索
  void set_spvec(size_t ncp, double m1, double sd1, double m2, double sd2);
  // 診断率配列をセットする
public:
  diag() : v_pi(BLANKf), v_uporder(true), covec(Nco, BLANKf), snvec(Nco, BLANKf), spvec(Nco, BLANKf), acvec(Nco, BLANKf), ppvvec(Nco, BLANKf), npvvec(Nco, BLANKf), v_d2(BLANKf), v_ppr(BLANKf), v_aucn(BLANKf), v_resdiag(resnon), statvec(2)
  // デフォールトコンストラクタ
  {}
  explicit diag(double pi0) : v_pi(pi0), v_uporder(true), covec(Nco, BLANKf), snvec(Nco, BLANKf), spvec(Nco, BLANKf), acvec(Nco, BLANKf), ppvvec(Nco, BLANKf), npvvec(Nco, BLANKf), v_d2(BLANKf), v_ppr(BLANKf), v_aucn(BLANKf), v_resdiag(resnon), statvec(2)
  // コンストラクタ
  {}
  virtual ~diag()               // デストラクタ
  {}
  diag(const diag& diag0) : v_pi(diag0.v_pi), v_uporder(diag0.v_uporder), covec(diag0.covec), snvec(diag0.snvec), spvec(diag0.spvec), acvec(diag0.acvec), ppvvec(diag0.ppvvec), npvvec(diag0.npvvec), v_d2(diag0.v_d2), v_ppr(diag0.v_ppr), v_aucn(diag0.v_aucn), v_resdiag(diag0.v_resdiag), statvec(diag0.statvec)
  // コピーコンストラクタ
  {}
  const diag& operator=(const diag& diag0); // コピー代入演算子
  double pi() const             // 群1の事前確率を返す
  { return v_pi; }
  double pi(double pi0)         // 群1の事前確率を設定する
  { return v_pi = pi0; }
  bool uporder() const          // データの昇順/下順の区別を返す
  { return v_uporder; }
  double co(size_t idx = 0) const // 理論的境界値を返す
  { return (idx < Nco) ? covec[idx] : 0.0; }
  double sn(size_t idx = 0) const // 理論的感度を返す
  { return (idx < Nco) ? snvec[idx] : 0.0; }
  double sp(size_t idx = 0) const // 理論的特異度を返す
  { return (idx < Nco) ? spvec[idx] : 0.0; }
  double ac(size_t idx = 0) const // 理論的正診率を返す 
  { return (idx < Nco) ? acvec[idx] : 0.0; }
  double ppv(size_t idx = 0) const // 理論的PPVを返す 
  { return (idx < Nco) ? ppvvec[idx] : 0.0; }
  double npv(size_t idx = 0) const // 理論的NPVを返す 
  { return (idx < Nco) ? npvvec[idx] : 0.0; }
  double d2() const             // マハラノビスの汎距離を返す
  { return v_d2; }
  double ppr() const            // 誤判別確率を返す
  { return v_ppr; }
  double aucn() const           // 理論的AUCを返す
  { return v_aucn; }
  double fsn(double z) const    // 理論的感度を返す
  { return (uporder()) ? 1.0 - (double)(sndis(z)) : (double)(sndis(z)); }
  double fsp(double z) const    // 理論的特異度を返す
  { return (uporder()) ? (double)(sndis(z)) : 1.0 - (double)(sndis(z)); }
  double fac(double z1, double z2) const // 理論的正診率を返す
  { return pi() * fsn(z1) + (1.0 - pi()) * fsp(z2); }
  double fppv(double sn, double sp) const; // 感度と特異度から陽性予測値を返す
  double fnpv(double sn, double sp) const; // 感度と特異度から陰性予測値を返す
  restype resdiag() const                  // 計算結果の種類を返す
  { return v_resdiag; }
  statv<double>& dstat(size_t idx) // 診断項目の基礎統計量を返す
  { return (idx < 2) ? statvec[idx] : statvec[0]; }
  void cal_dpa(option& opt, ldata& datalst); // DP-plot分析
  void cal_wtest(double& dev, double& vd, double& t, int& df); // Welchの検定
};

const diag& diag::operator=(const diag& diag0) // コピー代入演算子
{
  if (&diag0 == this) return *this; // 自分自身への代入を防ぐ
  v_pi = diag0.v_pi;  v_uporder = diag0.v_uporder;
  for (size_t i = 0; i < Nco; ++i)
  {
    covec[i] = diag0.covec[i];  snvec[i] = diag0.snvec[i];
    spvec[i] = diag0.spvec[i];  acvec[i] = diag0.acvec[i];
    ppvvec[i] = diag0.ppvvec[i];  npvvec[i] = diag0.npvvec[i];
  }
  v_d2 = diag0.v_d2;  v_ppr = diag0.v_ppr;  v_aucn = diag0.v_aucn;
  v_resdiag = diag0.v_resdiag;
  statvec.resize(diag0.statvec.size());
  statvec[0] = diag0.statvec[0];  statvec[1] = diag0.statvec[1];
  return *this;
}

double diag::fppv(double sn, double sp) const // 感度と特異度から陽性予測値を返す
{
  if (sgn(pi()) <= 0) return BLANKf;
  if (sgn(sn) <= 0) return uporder() ? 1.0 : 0.0;
  double ppv = 1.0 / (1.0 + (1.0 - pi()) / pi() * (1.0 - sp) / sn);
  if (sgn(ppv) <= 0) return 0.0;
  return (ppv > 1.0) ? 1.0 : ppv;
}

double diag::fnpv(double sn, double sp) const // 感度と特異度から陰性予測値を返す
{
  if (sgn(1.0 - pi()) <= 0) return BLANKf;
  if (sgn(sp) <= 0) return uporder() ? 1.0 : 0.0;
  double npv = 1.0 / (1.0 + pi() / (1.0 - pi()) * (1.0 - sn) / sp);
  if (sgn(npv) <= 0) return 0.0;
  return (npv > 1.0) ? 1.0 : npv;
}

void diag::search_acc(option& opt, ldata& datalst) // 実測値正診率点の探索
{
  int n1 = statvec[0].n(), n2 = statvec[1].n();
  int n = n1 + n2, tp = n1, dfn = 0, tn = 0, dtn = 0;
  double f = 0.0, fpre = BLANKf;
  int idx0 = 7, didx = 1;
  if (uporder() == false)        // 境界値が上限の時
  { idx0 = 8;  didx = -1; }
  int idx = idx0;
  for (auto& cdatai : datalst)
  {
    f = cdatai.fdata(1);
    if (f != fpre)
    {
      tp -= dfn;  tn += dtn;
      double sn = (double)tp / (double)n1;
      double sp = (double)tn / (double)n2;
      double ac = (double)(tp + tn) / (double)n;
      if (ac >= opt.ac())
      {
        covec[idx] = f;
        snvec[idx] = sn;  spvec[idx] = sp;  acvec[idx] = ac;
        ppvvec[idx] = (n2 - tn + tp > 0) ? (double)tp / (double)(n2 - tn + tp) : 0.0;
        npvvec[idx] = (tn + n1 - tp > 0) ? (double)tn / (double)(tn + n1 - tp) : 0.0;
        if (idx == idx0) idx += didx;
      }
      fpre = f;  dfn = dtn = 0;
    }
    if ((int)(cdatai.fdata(0)) == 0) ++dfn;
    else ++dtn;
  }
}

void diag::cal_uda(option& opt)  // 単変量線形判別分析
{
  int n1 = statvec[0].n(), n2 = statvec[1].n();
  double m1 = statvec[0].mean(), m2 = statvec[1].mean();
  double v1 = statvec[0].vr(), v2 = statvec[1].vr();
  double sd1 = statvec[0].sd(), sd2 = statvec[1].sd();
  if (v1 + v2 < FDeps) return;  // データが定数の時
  covec[0] = (m1 + m2) * 0.5;   // 判別分析の境界値
  if (sgndiff(m1, m2) == 0 || (sgndiff(pi(), 0.5) == 0 && sgndiff(sd1, sd2) == 0))
  { covec[1] = covec[2] = covec[3] = covec[0]; } // 全て判別分析の境界値にする
  else
  {
    covec[1] = cal_cospmin(m1, sd1, m2, sd2);
    covec[2] = cal_coacmax(m1, sd1, m2, sd2);
    covec[3] = cal_copnmin(m1, sd1, m2, sd2);
    if (opt.ac() != BLANKf)
    {
      double acmax = (covec[2] == BLANKf) ? 0.0 : fac((covec[2] - m1) / sd1, (covec[2] - m2) / sd2);
      if (opt.ac() <= acmax)
      {
        vecdbl covect(2);
        if (search_coacc(covect, opt.ac(), m1, sd1, m2, sd2) > 0)
        { covec[4] = covect[0];  covec[5] = covect[1]; }
      }
    }
  }
  v_d2 = 2.0 * (m1 - m2) * (m1 - m2) / (v1 + v2);
  double z = sqrt(v_d2 * 0.5);
  v_aucn = (double)(sndis(z));   // 理論AUC
  z = sqrt(v_d2) * 0.5;
  v_ppr = 1.0 - (double)(sndis(z)); // 誤判別確率
  if ((m1 < m2 && uporder()) || (m1 > m2 && uporder() == false)) v_ppr = 1.0 - v_ppr;
  covec[6] = opt.co();          // オプションで指定された境界値
//  set_spvec(Nco - 2, m1, sd1, m2, sd2); // 感度・特異度配列を設定する(最後の2つは設定しない)
  set_spvec(7, m1, sd1, m2, sd2); // オプションで設定された境界値まで設定する
}

double diag::cal_cospmin(double m1, double sd1, double m2, double sd2)
// 理論的|SN-SP|≒0の時の境界値を返す
{
  if (sd1 < FDeps || sd2 < FDeps) return (m1 + m2) * 0.5;
  return (m1 * sd2 + m2 * sd1) / (sd1 + sd2);
}

double diag::cal_coacmax(double m1, double sd1, double m2, double sd2)
// 理論的正診率最大値の時の境界値を計算する
{
  double v1 = sd1 * sd1, v2 = sd2 * sd2;
  double vd = v2 - v1;
  if (sgndiff(v2, v1) == 0)     // 2群の分散が等しい時
  {
    if (sgndiff(m1, m2) == 0) return (m1 + m2) * 0.5;
    return (m1 + m2) * 0.5 + v1 * log((1.0 - pi()) / pi()) / (m1 - m2);
  }
  if (sd1 < FDeps || sd2 < FDeps) return BLANKf; // 計算不可能
  double md = m1 * v2 - m2 * v1;
  double tmp = md * md - vd * (m1 * m1 * v2 - m2 * m2 * v1 + 2.0 * v1 * v2 * log((1.0 - pi()) / pi() * sd1 / sd2));
  double com = (sgn(tmp) >= 0 && sgn(vd) != 0) ? (md - sqrt(tmp)) / vd : BLANKf;
  double cop = (sgn(tmp) >= 0 && sgn(vd) != 0) ? (md + sqrt(tmp)) / vd : BLANKf;
  return (uporder()) ? com : cop; // 理論値を返す
}

double diag::cal_copnmin(double m1, double sd1, double m2, double sd2)
// 理論的|PPV-NPV|≒0の時の境界値を計算する
{
  if (sd1 < FDeps || sd2 < FDeps) return BLANKf; // 計算不可能
  constexpr int Rmax = 200;                      // 反復回数最大値
  constexpr double SDx = 5.0;   // 範囲を決めるためのSDの倍率
  double co = (m1 + m2) * 0.5;
  double z1 = (co - m1) / sd1, z2 = (co - m2) / sd2;
  double sn = fsn(z1), sp = fsp(z2);
  double pr = fppv(sn, sp) - fnpv(sn, sp);
  const double Delta = fabs(co) * 0.00001; // 刻み幅
  double comin = m2 - sd1 * SDx, comax = m1 + sd2 * SDx;
  if (m1 < m2)                  // 群1の平均値の方が小さい時
  { comin = m1 - sd1 * SDx;  comax = m2 + sd2 * SDx; }
  int nr = 0;
  while (fabs(pr) >= Rerror && nr < Rmax)
  {
    ++nr;
    double dco = Delta * co;
    double dco1 = dco / sd1, dco2 = dco / sd2;
    double snd = fsn(z1 + dco1), spd = fsp(z2 + dco2);
    double prd = fppv(snd, spd) - fnpv(snd, spd);
    double slope = (prd - pr) / dco;
    if (fabs(slope) < FDeps) break; // 収束
    co -= pr / slope;
    if (co < comin || comax < co) // 範囲外
    {
      co = BLANKf;
      break;
    }
    z1 = (co - m1) / sd1;  z2 = (co - m2) / sd2;
    sn = fsn(z1);  sp = fsp(z2);
    pr = fppv(sn, sp) - fnpv(sn, sp);
  }
  if (co == BLANKf) return co;  // 計算できない時
  return co;                    // 境界値を返す
}

int diag::search_coacc(vecdbl& covect, double acc, double m1, double sd1, double m2, double sd2)
// 理論的正診率点の探索
{
  constexpr double SDx = 4.0;         // SDの倍数
  if (sd1 < FDeps || sd2 < FDeps) return 0; // データが定数の時は計算できず
  double comin = m2 - sd2 * SDx, comax = m1 + sd1 * SDx;
  if (m1 < m2)                  // 群1の平均値の方が小さい時
  { comin = m1 - sd1 * SDx;  comax = m2 + sd2 * SDx; }
  constexpr uint RDdef = 1000;  // 分割数
  const double Dco = (comax - comin) / (double)RDdef;
  covect[0] = covect[1] = BLANKf;
  int idx = 0;
  double co = comin;
  for (uint i = 0; i <= RDdef; ++i, co += Dco)
  {
    double ac = fac((co - m1) / sd1, (co - m2) / sd2);
    if (ac >= acc)
    {
      covect[idx] = co;
      if (idx == 0) ++idx;
    }
  }
  int nco = 0;
  if (covect[0] == BLANKf) covect[0] = comin;
  else nco = 1;
  if (covect[1] == BLANKf) covect[1] = comax;
  else ++nco;
  return nco;                   // 境界値数を返す
}

void diag::set_spvec(size_t ncp, double m1, double sd1, double m2, double sd2)
// 診断率配列をセットする
{
  for (size_t i = 0; i < ncp; ++i)
  {
    if (covec[i] != BLANKf)
    {
      if (fabs(covec[i]) < FDeps) covec[i] = 0.0;
      double sn = 0.0, sp = 0.0, z1 = 0.0, z2 = 0.0;
      if (sd1 > FDeps)
      {
        z1 = (covec[i] - m1) / sd1;
        sn = fsn(z1);
      }
      else
      {
        z1 = 0.0;
        if (uporder()) sn = (m1 < covec[i]) ? 0.0 : 1.0;
        else sn = (covec[i] < m1) ? 0.0 : 1.0;
      }
      if (sd2 > FDeps)
      {
        z2 = (covec[i] - m2) / sd2;
        sp = fsp(z2);
      }
      else
      {
        z2 = 0.0;
        if (uporder()) sp = (covec[i] < m2) ? 0.0 : 1.0;
        else sp = (m2 < covec[i]) ? 0.0 : 1.0;
      }
      snvec[i] = sn;  spvec[i] = sp;
      acvec[i] = pi() * sn + (1.0 - pi()) * sp;
      ppvvec[i] = fppv(sn, sp);  npvvec[i] = fnpv(sn, sp);
    }
  }
}

void diag::cal_dpa(option& opt, ldata& datalst) // DP-plot分析
{
  for (auto& cdatai : datalst) statvec[(size_t)(cdatai.fdata(0))].sigma(cdatai.fdata(1)); // 群ごとに統計量を集計する
  int n1 = statvec[0].n(), n2 = statvec[1].n();
  if (n1 * n2 < 1)              // 計算できず
  {
    v_resdiag = reserr;
    return;
  }
  int nt = n1 + n2;
  if (opt.lim() == limlow)      // 境界値が下限の時
  { v_uporder = true; }
  else if (opt.lim() == limup)  // 境界値が上限の時
  { v_uporder = false; }
  else                          // 自動判定
  { v_uporder = (statvec[0].mean() >= statvec[1].mean()) ? true : false; }
  if (v_uporder == false) datalst.reverse(); // データリストを逆順にする
  if (v_pi == BLANKf) v_pi = (double)n1 / (double)nt; // v_piを現在のデータの割合で設定する
  if (opt.ac() != BLANKf) search_acc(opt, datalst);
  cal_uda(opt);
  v_resdiag = resok;             // DP-plot分析成功
}

void diag::cal_wtest(double& dev, double& vd, double& t, int& df)
// Welchの検定結果を返す
{
  dev = fabs(statvec[0].mean() - statvec[1].mean());
  int n1 = statvec[0].n(), n2 = statvec[1].n();
  double v1 = statvec[0].vr() / (double)n1;
  double v2 = statvec[1].vr() / (double)n2;
  vd = v1 + v2;
  df = (vd > 0.0) ? (int)(vd * vd / (v1 * v1 / (n1 - 1) + v2 * v2 / (n2 -1))) : n1 + n2  - 2;
  t = (vd > 0.0) ? dev / sqrt(vd) : 0.0;
}

class dpplotclass               // dpplotメインクラス
{
private:
  string fname;                 // ファイル名
  void put_resid(ostream& fout, idd& id) const;
  // 分析結果のタイトルとデータ識別情報を出力する
  void put_result(ostream& fout, option& opt, idd& id, ldata& datalst, diag& dg) const;
  // 分析結果を出力する
  void put_stat(ostream& fout, option& opt, diag& dg) const;
  // 基礎統計量を出力する
  void put_dpa(ostream& fout, option& opt, diag& dg) const;
  // DP-plot分析の結果を出力する
  void put_acc(ostream& fout, option& opt, const diag& dg) const;
  // 指定した正診率になる時の診断率を出力する
  void put_uda(ostream& fout, option& opt, diag& dg) const;
  // 単変量線型判別分析の結果とWelchの検定結果を出力する
  void plot_ngp(option& opt, idd& id, ldata& datalst, diag& dg) const;
  // Ngraph用データを出力する
  void plot_ngpid(ofstream& fout, idd& id) const;
  // Ngraph用データのタイトルとデータ識別情報を出力する
  void plot_ac(ofstream& fout, option& opt, ldata& datalst, diag& dg) const;
  // Ngraph用正診率曲線を出力する
  void plot_acn(ofstream& fout, option& opt, diag& dg) const;
  // Ngraph用理論的正診率曲線を出力する
  double set_xrange(double& xmin, double& xmax, option& opt, diag& dg) const;
  // X軸の最小値と最大値を設定して範囲を返す
  void plot_log(ofstream& fout, option& opt, diag& dg) const;
  // Ngraph用ロジスティック曲線を出力する
  void plot_logci(ofstream& fout, option& opt, diag& dg) const;
  // Ngraph用ロジスティック曲線の信頼区間を出力する
  void plot_pnfn(ofstream& fout, option& opt, diag& dg) const;
  // Ngraph用理論的陽性率/陰性率曲線を出力する 
  void plot_pv(ofstream& fout, ldata& datalst, diag& dg) const;
  // Ngraph用的PPV/NPV曲線を出力する
  void plot_pvn(ofstream& fout, option& opt, diag& dg) const;
  // Ngraph用理論的PPV/NPV曲線を出力する
  void plot_pvnci(ofstream& fout, option& opt, diag& dg) const;
  // Ngraph用理論的PPV/NPV曲線の信頼区間を出力する
  void plot_snp(ofstream& fout, option& opt, ldata& datalst, diag& dg) const;
  // Ngraph用感度/特異度曲線を出力する
  void plot_snpn(ofstream& fout, option& opt, diag& dg) const;
  // Ngraph用理論的感度/特異度曲線を出力する
public:
  dpplotclass() : fname()       // コンストラクタ
  {}
  ~dpplotclass()                // デストラクタ
  {}
  int get(const char* c_fname, option& opt, idd& id, ldata& datalst);
  // データを読み込み累積する
  void put(option& opt, idd& id, ldata& datalst) const; // 集計結果を出力する
};

void dpplotclass::put_resid(ostream& fout, idd& id) const
// 分析結果のタイトルとデータ識別情報を出力する
{
  puttitle(fout, Title);
  fout << "\n";
  id.put_title(fout);
  id.put_stra(fout);
  fout << "群項目  :" << id.item(0) << "\n";
  fout << "集計項目:" << id.item(1) << "\n";
}

void dpplotclass::put_result(ostream& fout, option& opt, idd& id, ldata& datalst, diag& dg) const
// 分析結果を出力する
{
  constexpr int LENline = 74;   // 区切り線の長さ
  put_resid(fout, id);
  bool putres = false;
  if (opt.res("stat"))
  {
    fout << "\n";
    putline(fout, LENline);
    put_stat(fout, opt, dg);
    putres = true;
  }
  if (opt.res("dpa"))
  {
    if (putres == false) fout << "\n";
    putline(fout, LENline);
    put_dpa(fout, opt, dg);
    putres = true;
  }
  if (opt.res("ac"))
  {
    if (putres == false) fout << "\n";
    putline(fout, LENline);
    put_acc(fout, opt, dg);
    putres = true;
  }
  if (opt.res("uda"))
  {
    if (putres == false) fout << "\n";
    putline(fout, LENline);
    put_uda(fout, opt, dg);
    putres = true;
  }
  if (putres) putline(fout, LENline);
  fout << "\n" << endl;
}

void dpplotclass::put_stat(ostream& fout, option& opt, diag& dg) const
// 基礎統計量を出力する
{
  fout << "群1(疾患群):コード=";
  if ((opt.g(0)).empty()) fout << opt.g(1) << " 以外\n";
  else fout << opt.g(0) << "\n";
  fout << "  ";
  putmsd(fout, (dg.dstat(0)).n(), (dg.dstat(0)).mean(), (dg.dstat(0)).vr());
  fout << "群2(正常群):コード=";
  if ((opt.g(1)).empty()) fout << opt.g(0) << " 以外\n";
  else fout << opt.g(1) << "\n";
  fout << "  ";
  putmsd(fout, (dg.dstat(1)).n(), (dg.dstat(1)).mean(), (dg.dstat(1)).vr());
}

void dpplotclass::put_dpa(ostream& fout, option& opt, diag& dg) const
// DP-plot分析の結果を出力する
{
  fout << "SS-point(理論的感度≒理論的特異度の点)=";
  if (dg.co(1) == BLANKf)
  { fout << "なし\n"; }
  else
  {
    string tstr("以上");
    if (dg.uporder() == false) tstr = "以下";
    fout << dg.co(1) << tstr << "を群1と診断した場合\n";
    fout << "  理論的感度=" << dg.sn(1) << "  理論的特異度=" << dg.sp(1) << "  理論的正診率=" << dg.ac(1) << "\n";
    fout << "  理論的陽性予測値=" << dg.ppv(1) << "  理論的陰性予測値=" << dg.npv(1) << "\n";
  }
  fout << "AC-point(理論的正診率最大点)=";
  if (dg.co(2) == BLANKf)
  { fout << "なし\n"; }
  else
  {
    fout << dg.co(2) << "(群1の事前確率=" << dg.pi() << ")\n";
    fout << "  理論的感度=" << dg.sn(2) << "  理論的特異度=" << dg.sp(2) << "  理論的正診率=" << dg.ac(2) << "\n";
    fout << "  理論的陽性予測値=" << dg.ppv(2) << "  理論的陰性予測値=" << dg.npv(2) << "\n";
  }
  fout << "PN-point(理論的PPV≒理論的NPVの点)=";
  if (dg.co(3) == BLANKf)
  { fout << "なし\n"; }
  else
  {
    fout << dg.co(3) << "(群1の事前確率=" << dg.pi() << ")\n";
    fout << "  理論的感度=" << dg.sn(3) << "  理論的特異度=" << dg.sp(3) << "  理論的正診率=" << dg.ac(3) << "\n";
    fout << "  理論的陽性予測値=" << dg.ppv(3) << "  理論的陰性予測値=" << dg.npv(3) << "\n";
  }
  if (opt.co() == BLANKf) return;
  fout << "指定境界値=" << dg.co(6) << "(群1の事前確率=" << dg.pi() << ")\n";
  fout << "  理論的感度=" << dg.sn(6) << "  理論的特異度=" << dg.sp(6) << "  理論的正診率=" << dg.ac(6) << "\n";
  fout << "  理論的陽性予測値=" << dg.ppv(6) << "  理論的陰性予測値=" << dg.npv(6) << endl;
}

void dpplotclass::put_acc(ostream& fout, option& opt, const diag& dg) const
// 指定した正診率になる時の診断率を出力する
{
  fout << "理論的正診率が" << opt.ac() << "以上になる最小値点=";
  if (dg.co(4) == BLANKf)
  { fout << "なし\n"; }
  else
  {
    fout << dg.co(4) << "\n";
    fout << "  理論的感度=" << dg.sn(4) << "  理論的特異度=" << dg.sp(4) << "  理論的正診率=" << dg.ac(4) << "\n";
    fout << "  理論的陽性予測値=" << dg.ppv(4) << "  理論的陰性予測値=" << dg.npv(4) << "\n";
  }
  fout << "理論的正診率が" << opt.ac() << "以上になる最大値点=";
  if (dg.co(5) == BLANKf)
  { fout << "なし\n"; }
  else
  {
    fout << dg.co(5) << "\n";
    fout << "  理論的感度=" << dg.sn(5) << "  理論的特異度=" << dg.sp(5) << "  理論的正診率=" << dg.ac(5) << "\n";
    fout << "  理論的陽性予測値=" << dg.ppv(5) << "  理論的陰性予測値=" << dg.npv(5) << "\n";
  }
  fout << "正診率が" << opt.ac() << "以上になる最小値点=";
  if (dg.co(7) == BLANKf)
  { fout << "なし\n"; }
  else
  {
    fout << dg.co(7) << "\n";
    fout << "  感度=" << dg.sn(7) << "  特異度=" << dg.sp(7) << "  正診率=" << dg.ac(7) << "\n";
    fout << "  陽性予測値=" << dg.ppv(7) << "  陰性予測値=" << dg.npv(7) << "\n";
  }
  fout << "正診率が" << opt.ac() << "以上になる最大値点=";
  if (dg.co(8) == BLANKf)
  { fout << "なし" << endl; }
  else
  {
    fout << dg.co(8) << "\n";
    fout << "  感度=" << dg.sn(8) << "  特異度=" << dg.sp(8) << "  正診率=" << dg.ac(7) << "\n";
    fout << "  陽性予測値=" << dg.ppv(8) << "  陰性予測値=" << dg.npv(8) << endl;
  }
}

void dpplotclass::put_uda(ostream& fout, option& opt, diag& dg) const
// 単変量線型判別分析の結果とWelchの検定結果を出力する
{
  string tstr("以上");
  if (dg.uporder() == false) tstr = "以下";
  fout << "単変量線型判別分析:境界値=" << dg.co(0) << tstr << "を群1と診断した場合\n";
  fout << "  理論的感度=" << dg.sn(0) << "  理論的特異度=" << dg.sp(0) << "  理論的正診率=" << dg.ac(0) << "\n";
  fout << "  理論的陽性予測値=" << dg.ppv(0) << "  理論的陰性予測値=" << dg.npv(0) << "\n";
  double n1 = (double)((dg.dstat(0)).n()), n2 = (double)((dg.dstat(1)).n());
  double f = dg.d2() * n1 * n2 / (n1 + n2);
  double r2 = f / (f + n1 + n2 - 2.0);
  fout << "判別効率(Mahalanobisの平方距離):D2乗=" << dg.d2() << "  相関比:η=" << sqrt(r2) << "\n";
  fout << "  判別確率:Φ(D/2)=Φ(" << sqrt(dg.d2()) * 0.5 << ")=" << 1.0 - dg.ppr() << "(誤判別確率:1-p=" << dg.ppr() << ")\n";
  fout << "  理論的ROC曲線のAUC=Φ(D/√2)=Φ(" << sqrt(dg.d2() * 0.5) << ")=" << dg.aucn() << "\n";
  double dev = 0.0, vd = 0.0, t = 0.0;
  int df = 0;
  dg.cal_wtest(dev, vd, t, df);
  ldbl pr = ttop(t, df);
  fout << "Welchの検定:t=" << t << "  自由度=" << df << "  有意確率:p=" << pr << pind(pr) << "\n";
  double ci = calmcci(mcfis, opt.cc(), 2, 1, df, 1, vd);
  fout << "  平均値の差の" << opt.cc() * 100.0 << "%CI=";
  fout << dev << "±" << ci << "(" << dev - ci << "-" << dev + ci << ")\n";
}

void dpplotclass::plot_ngp(option& opt, idd& id, ldata& datalst, diag& dg) const
// Ngraph用データを出力する
{
  ofstream fout;
  if (ofopen(fout, (opt.ngpname()).c_str()) == false) return; // ファイルオープン失敗
  bool delml = false;
  if (opt.ngp("dpp") || opt.ngp("snp")) // 感度・特異度曲線を出力する
  {
    plot_ngpid(fout, id);
    plot_snp(fout, opt, datalst, dg);
    fout << NGPdelm << "\n";
    plot_snpn(fout, opt, dg);
    delml = true;
  }
  if (opt.ngp("dpp") || opt.ngp("ac")) // 正診率曲線を出力する
  {
    if (delml)
    { fout << NGPdelm << "\n"; }
    else
    {
      fout << "-s1 -dline -w80 -zy0,1,0.2\n";
      delml = true;
    }
    plot_ngpid(fout, id);
    if (opt.ngp("ac"))
    {
      plot_ac(fout, opt, datalst, dg);
      fout << NGPdelm << "\n";
    }
    plot_acn(fout, opt, dg);
  }
  if (opt.ngp("pv"))            // PPV/NPV曲線を出力する
  {
    if (delml)
    { fout << NGPdelm << "\n"; }
    else
    {
      fout << "-s1 -dline -w80 -zy0,1,0.2\n";
      delml = true;
    }
    plot_ngpid(fout, id);
    plot_pv(fout, datalst, dg);
    fout << NGPdelm << "\n";
    plot_pvn(fout, opt, dg);
    fout << NGPdelm << "\n";
    plot_pvnci(fout, opt, dg);
  }
  if (opt.ngp("pnf"))           // 陽性率／陰性率曲線を出力する
  {
    if (delml)
    { fout << NGPdelm << "\n"; }
    else
    {
      fout << "-s1 -dline -w80 -zy0,1,0.2\n";
      delml = true;
    }
    plot_ngpid(fout, id);
    plot_pnfn(fout, opt, dg);
  }
  if (opt.ngp("log"))           // ロジスティック曲線を出力する
  {
    if (delml) fout << NGPdelm << "\n";
    else fout << "-s1 -dline -w80 -zy0,1,0.2\n";
    plot_ngpid(fout, id);
    plot_log(fout, opt, dg);
    fout << NGPdelm << "\n";
    plot_logci(fout, opt, dg);
  }
}

void dpplotclass::plot_ngpid(ofstream& fout, idd& id) const
// Ngraph用タイトルとデータ識別情報を出力する
{
  fout << "%";  puttitle(fout, Title);
  fout << "%\n";
  id.put_title(fout, formngp);
  id.put_stra(fout, formngp);
  fout << "%群項目  :" << id.item(0) << "\n";
  fout << "%集計項目:" << id.item(1) << "\n";
}

void dpplotclass::plot_ac(ofstream& fout, option& opt, ldata& datalst, diag& dg) const
// Ngraph用正診率曲線を出力する
{
  int n1 = (dg.dstat(0)).n(), n2 = (dg.dstat(1)).n();
  int tp = n1, tn = 0;
  double n = (double)(n1 + n2);
  double pi1 = (double)n1 / (double)n;
  fout << "%正診率曲線:群1の事前確率=" << pi1 << "\n";
  double f = 0.0, ac = 0.0;
  for (auto& cdatai : datalst)
  {
    f = cdatai.fdata(1);
    ac = (double)(tp + tn) / n;
    fout << f << " " << ac << "\n";
    if ((int)(cdatai.fdata(0)) == 0) --tp;
    else ++tn;
  }
  ac = (double)n2 / n;
  fout << f << " " << ac << "\n";
  fout << NGPdelm2 << "\n";
  pi1 = (0.0 < opt.pi() && opt.pi() < 1.0) ? opt.pi() : 0.5;
  fout << "%正診率曲線:群1の事前確率=" << pi1 << "\n"; // 指定された事前確率
  tp = n1;  tn = 0;
  for (auto& cdatai : datalst)
  {
    f = cdatai.fdata(1);
    double sn = (double)tp / (double)n1;
    double sp = (double)tn / (double)n2;
    ac = pi1 * sn + (1.0 - pi1) * sp;
    fout << f << " " << ac << "\n";
    if ((int)(cdatai.fdata(0)) == 0) --tp;
    else ++tn;
  }
  ac = 1.0 - pi1;
  fout << f << " " << ac << endl;
}

void dpplotclass::plot_acn(ofstream& fout, option& opt, diag& dg) const
// Ngraph用理論的正診率曲線を出力する
{
  int n1 = (dg.dstat(0)).n(), n2 = (dg.dstat(1)).n();
  double m1 = (dg.dstat(0)).mean(), sd1 = (dg.dstat(0)).sd();
  double m2 = (dg.dstat(1)).mean(), sd2 = (dg.dstat(1)).sd();
  double xmin = 0.0, xmax = 0.0;
  double xrange = set_xrange(xmin, xmax, opt, dg);
  double dx = xrange / (double)RDdef;
  double pi0 = dg.pi();         // 現在の事前確率を保存
  double pi1 = (double)n1 / (double)(n1 + n2);
  dg.pi(pi1);
  fout << "%理論的正診率曲線:群1の事前確率=" << pi1 << "\n";
  double x = xmin;
  for (uint i = 0; i <= RDdef; ++i, x += dx)
  { fout << x << " " << dg.fac((x - m1) / sd1, (x - m2) / sd2) << "\n"; }
  fout << NGPdelm2 << "\n";
  pi1 = (0.0 < opt.pi() && opt.pi() < 1.0) ? opt.pi() : 0.5;
  dg.pi(pi1);
  fout << "%理論的正診率曲線:群1の事前確率=" << pi1 << "\n"; // 指定された事前確率
  x = xmin;
  for (uint i = 0; i <= RDdef; ++i, x += dx)
  { fout << x << " " << dg.fac((x - m1) / sd1, (x - m2) / sd2) << "\n"; }
  dg.pi(pi0);                   // 事前確率を元に戻す
}

double dpplotclass::set_xrange(double& xmin, double& xmax, option& opt, diag& dg) const
// X軸の最小値と最大値を設定して範囲を返す
{
  if (opt.xmin() != BLANKf && opt.xmax() != BLANKf) // X軸の範囲指定あり
  { xmin = opt.xmin();  xmax = opt.xmax(); }
  else
  {
    constexpr double SDx = 4.0;   // SDの倍数
    int n1 = (dg.dstat(0)).n(), n2 = (dg.dstat(1)).n();
    double m1 = (dg.dstat(0)).mean(), sd1 = (dg.dstat(0)).sd();
    double m2 = (dg.dstat(1)).mean(), sd2 = (dg.dstat(1)).sd();
    xmin = m2 - sd2 * SDx;  xmax = m1 + sd1 * SDx;
    if (m1 < m2)                  // 群1の平均値の方が小さい時
    { xmin = m1 - sd1 * SDx;  xmax = m2 + sd2 * SDx; }
  }
  return xmax - xmin;           // 範囲を返す
}

void dpplotclass::plot_log(ofstream& fout, option& opt, diag& dg) const
// Ngraph用ロジスティック曲線を出力する
{
  int n1 = (dg.dstat(0)).n(), n2 = (dg.dstat(1)).n();
  double m1 = (dg.dstat(0)).mean(), sd1 = (dg.dstat(0)).sd();
  double m2 = (dg.dstat(1)).mean(), sd2 = (dg.dstat(1)).sd();
  double xmin = 0.0, xmax = 0.0;
  double xrange = set_xrange(xmin, xmax, opt, dg);
  double dx = xrange / (double)RDdef;
  double orpi = (1.0 - dg.pi()) / dg.pi();
  fout << "%ロジスティック曲線:群1の事前確率=" << dg.pi() << "\n";
  double x = xmin;
  for (uint i = 0; i <= RDdef; ++i, x += dx)
  {
    double p1 = snden((x - m1) / sd1), p2 = snden((x - m2) / sd2);
    double pr = 1.0 / (1.0 + orpi * p2 / p1);
    if (dg.uporder() == false) pr = 1.0 - pr;
    fout << x << " " << pr << "\n";
  }
}

void dpplotclass::plot_logci(ofstream& fout, option& opt, diag& dg) const
// Ngraph用ロジスティック曲線の信頼区間を出力する
{
  int n1 = (dg.dstat(0)).n(), n2 = (dg.dstat(1)).n();
  double m1 = (dg.dstat(0)).mean(), sd1 = (dg.dstat(0)).sd();
  double m2 = (dg.dstat(1)).mean(), sd2 = (dg.dstat(1)).sd();
  double xmin = 0.0, xmax = 0.0;
  double xrange = set_xrange(xmin, xmax, opt, dg);
  double dx = xrange / (double)RDdef;
  double orpi = (1.0 - dg.pi()) / dg.pi();
  double t = ptoz(1.0 - opt.cc());
  fout << "%ロジスティック曲線の" << opt.cc() * 100.0 << "%CI下限:群1の事前確率=" << dg.pi() << "\n";
  double x = xmin;
  for (uint i = 0; i <= RDdef; ++i, x += dx)
  {
    double p1 = snden((x - m1) / sd1), p2 = snden((x - m2) / sd2);
    double p = 1.0 / (1.0 + orpi * p2 / p1);
    double pl = p - t * sqrt(p * (1.0 - p) / (double)(n1 + n2));
    if (pl < 0.0) pl = 0.0;
    if (dg.uporder() == false) pl = 1.0 - pl;
    fout << x << " " << pl << "\n";
  }
  fout << NGPdelm2 << "\n";
  fout << "%ロジスティック曲線の" << opt.cc() * 100.0 << "%CI上限:群1の事前確率=" << dg.pi() << "\n";
  x = xmin;
  for (uint i = 0; i <= RDdef; ++i, x += dx)
  {
    double p1 = snden((x - m1) / sd1), p2 = snden((x - m2) / sd2);
    double p = 1.0 / (1.0 + orpi * p2 / p1);
    double pu = p + t * sqrt(p * (1.0 - p) / (double)(n1 + n2));
    if (pu > 1.0) pu = 1.0;
    if (dg.uporder() == false) pu = 1.0 - pu;
    fout << x << " " << pu << "\n";
  }
}

void dpplotclass::plot_pnfn(ofstream& fout, option& opt, diag& dg) const
// Ngraph用理論的陽性率/陰性率曲線を出力する
{
  int n1 = (dg.dstat(0)).n(), n2 = (dg.dstat(1)).n();
  double m1 = (dg.dstat(0)).mean(), sd1 = (dg.dstat(0)).sd();
  double m2 = (dg.dstat(1)).mean(), sd2 = (dg.dstat(1)).sd();
  double xmin = 0.0, xmax = 0.0;
  double xrange = set_xrange(xmin, xmax, opt, dg);
  double dx = xrange / (double)RDdef;
  fout << "%理論的陽性率曲線:群1の事前確率=" << dg.pi() << "\n";
  double x = xmin;
  for (uint i = 0; i <= RDdef; ++i, x += dx)
  {
    double sn = dg.fsn((x - m1) / sd1);
    double sp = dg.fsp((x - m2) / sd2);
    double pf = dg.pi() * sn + (1.0 - dg.pi()) * (1.0 - sp);
    fout << x << " " << pf << "\n";
  }
  fout << NGPdelm2 << "\n";
  fout << "%理論的陰性率曲線:群2の事前確率=" << 1.0 - dg.pi() << "\n";
  x = xmin;
  for (uint i = 0; i <= RDdef; ++i, x += dx)
  {
    double sn = dg.fsn((x - m1) / sd1);
    double sp = dg.fsp((x - m2) / sd2);
    double nf = dg.pi() * (1.0 - sn) + (1.0 - dg.pi()) * sp;
    fout << x << " " << nf << "\n";
  }
  fout << NGPdelm2 << "\n";
  fout << "%理論的疾患率曲線:群1の事前確率=" << dg.pi() << "\n";
  x = xmin;
  for (uint i = 0; i <= RDdef; ++i, x += dx)
  {
    double sn = dg.fsn((x - m1) / sd1);
    fout << x << " " << dg.pi() * sn << "\n";
  }
  fout << NGPdelm2 << "\n";
  fout << "%理論的正常率曲線:群2の事前確率=" << 1.0 - dg.pi() << "\n";
  x = xmin;
  for (uint i = 0; i <= RDdef; ++i, x += dx)
  {
    double sp = dg.fsp((x - m2) / sd2);
    fout << x << " " << (1.0 - dg.pi()) * sp << "\n";
  }
}

void dpplotclass::plot_pv(ofstream& fout, ldata& datalst, diag& dg) const
// Ngraph用PPV/NPV曲線を出力する
{
  int n1 = (dg.dstat(0)).n(), n2 = (dg.dstat(1)).n();
  int n = n1 + n2, tp = n1, tn = 0;
  double f = 0.0, ppv = 0.0;
  fout << "%PPV曲線:群1の事前確率=" << dg.pi() << "\n";
  for (auto& cdatai : datalst)
  {
    f = cdatai.fdata(1);
    double sn = (double)tp / (double)n1;
    double sp = (double)tn / (double)n2;
    ppv = dg.fppv(sn, sp);
    fout << f << " " << ppv << "\n";
    if ((int)(cdatai.fdata(0)) == 0) --tp;
    else ++tn;
  }
  tp = n1;  tn = 0;
  double npv = 0.0;
  fout << NGPdelm2 << "\n";
  fout << "%NPV曲線:群2の事前確率=" << 1.0 - dg.pi() << "\n";
  for (auto& cdatai : datalst)
  {
    f = cdatai.fdata(1);
    double sn = (double)tp / (double)n1;
    double sp = (double)tn / (double)n2;
    npv = dg.fnpv(sn, sp);
    fout << f << " " << npv << "\n";
    if ((int)(cdatai.fdata(0)) == 0) --tp;
    else ++tn;
  }
}

void dpplotclass::plot_pvn(ofstream& fout, option& opt, diag& dg) const
// Ngraph用理論的PPV/NPV曲線を出力する
{
  int n1 = (dg.dstat(0)).n(), n2 = (dg.dstat(1)).n();
  double m1 = (dg.dstat(0)).mean(), sd1 = (dg.dstat(0)).sd();
  double m2 = (dg.dstat(1)).mean(), sd2 = (dg.dstat(1)).sd();
  double xmin = 0.0, xmax = 0.0;
  double xrange = set_xrange(xmin, xmax, opt, dg);
  double dx = xrange / (double)RDdef;
  fout << "%理論的PPV:群1の事前確率=" << dg.pi() << "\n";
  double x = xmin;
  for (uint i = 0; i <= RDdef; ++i, x += dx)
  {
    double sn = dg.fsn((x - m1) / sd1);
    double sp = dg.fsp((x - m2) / sd2);
    double ppv = dg.fppv(sn, sp);
    fout << x << " " << ppv << "\n";
  }
  fout << NGPdelm2 << "\n";
  fout << "%理論的NPV:群2の事前確率=" << 1.0 - dg.pi() << "\n";
  x = xmin;
  for (uint i = 0; i <= RDdef; ++i, x += dx)
  {
    double sn = dg.fsn((x - m1) / sd1);
    double sp = dg.fsp((x - m2) / sd2);
    double npv = dg.fnpv(sn, sp);
    fout << x << " " << npv << "\n";
  }
}

void dpplotclass::plot_pvnci(ofstream& fout, option& opt, diag& dg) const
// Ngraph用理論的PPV/NPV曲線の信頼区間を出力する
{
   int n1 = (dg.dstat(0)).n(), n2 = (dg.dstat(1)).n();
  double m1 = (dg.dstat(0)).mean(), sd1 = (dg.dstat(0)).sd();
  double m2 = (dg.dstat(1)).mean(), sd2 = (dg.dstat(1)).sd();
  double xmin = 0.0, xmax = 0.0;
  double xrange = set_xrange(xmin, xmax, opt, dg);
  double dx = xrange / (double)RDdef;
  double t = ptoz(1.0 - opt.cc());
  fout << "%理論的PPVの" << opt.cc() * 100.0 << "%CI下限:群1の事前確率=" << dg.pi() << "\n";
  double x = xmin;
  for (uint i = 0; i <= RDdef; ++i, x += dx)
  {
    double sn = dg.fsn((x - m1) / sd1);
    double sp = dg.fsp((x - m2) / sd2);
    double ppv = dg.fppv(sn, sp);
    double ppvl = ppv - t * sqrt(ppv * (1.0 - ppv) / (double)(n1 + n2));
    if (ppvl < 0.0) ppvl = 0.0;
    fout << x << " " << ppvl << "\n";
  }
  fout << NGPdelm2 << "\n";
  fout << "%理論的PPVの" << opt.cc() * 100.0 << "%CI上限:群1の事前確率=" << dg.pi() << "\n";
  x = xmin;
  for (uint i = 0; i <= RDdef; ++i, x += dx)
  {
    double sn = dg.fsn((x - m1) / sd1);
    double sp = dg.fsp((x - m2) / sd2);
    double ppv = dg.fppv(sn, sp);
    double ppvu = ppv + t * sqrt(ppv * (1.0 - ppv) / (double)(n1 + n2));
    if (ppvu > 1.0) ppvu = 1.0;
    fout << x << " " << ppvu << "\n";
  }
  fout << NGPdelm2 << "\n";
  fout << "%理論的NPVの" << opt.cc() * 100.0 << "%CI下限:群2の事前確率=" << 1.0 - dg.pi() << "\n";
  x = xmin;
  for (uint i = 0; i <= RDdef; ++i, x += dx)
  {
    double sn = dg.fsn((x - m1) / sd1);
    double sp = dg.fsp((x - m2) / sd2);
    double npv = dg.fnpv(sn, sp);
    double npvl = npv - t * sqrt(npv * (1.0 - npv) / (double)(n1 + n2));
    if (npvl < 0.0) npvl = 0.0;
    fout << x << " " << npvl << "\n";
  }
  fout << NGPdelm2 << "\n";
  fout << "%理論的NPVの" << opt.cc() * 100.0 << "%CI上限:群2の事前確率=" << 1.0 - dg.pi() << "\n";
  x = xmin;
  for (uint i = 0; i <= RDdef; ++i, x += dx)
  {
    double sn = dg.fsn((x - m1) / sd1);
    double sp = dg.fsp((x - m2) / sd2);
    double npv = dg.fnpv(sn, sp);
    double npvu = npv + t * sqrt(npv * (1.0 - npv) / (double)(n1 + n2));
    if (npvu > 1.0) npvu = 1.0;
    fout << x << " " << npvu << "\n";
  }
}

void dpplotclass::plot_snp(ofstream& fout, option& opt, ldata& datalst, diag& dg) const
// Ngraph用感度/特異度曲線を出力する
{
  int n1 = (dg.dstat(0)).n(), n2 = (dg.dstat(1)).n();
  int n = n1 + n2, tp = n1, tn = 0;
  fout << "%感度曲線 群1:コード=";
  if ((opt.g(0)).empty()) fout << opt.g(1) << " 以外\n";
  else fout << opt.g(0) << "\n";
  double f = 0.0;
  for (auto& cdatai : datalst)
  {
    if ((int)(cdatai.fdata(0)) == 0) // 群1
    {
      f = cdatai.fdata(1);
      double sn = (double)tp / (double)n1;
      fout << f << " " << sn << "\n";
      --tp;
    }
  }
  fout << "%" << f << " 0\n"; // 感度0の点は実測値ではないのでプロットしない
  fout << NGPdelm2 << "\n";
  fout << "%特異度曲線 群2:コード=";
  if ((opt.g(1)).empty()) fout << opt.g(0) << " 以外\n";
  else fout << opt.g(1) << "\n";
  tp = n1;  tn = 0;
  for (auto& cdatai : datalst)
  {
    if ((int)(cdatai.fdata(0)) == 1) // 群2
    {
      f = cdatai.fdata(1);
      double sp = (double)tn / (double)n2;
      fout << f << " " << sp << "\n"; 
      ++tn;
    }
  }
  fout << "%" << f << " 1\n"; // 特異度1の点は実測値ではないのでプロットしない
}

void dpplotclass::plot_snpn(ofstream& fout, option& opt, diag& dg) const
// Ngraph用理論的感度／特異度曲線を出力する
{
  int n1 = (dg.dstat(0)).n(), n2 = (dg.dstat(1)).n();
  double m1 = (dg.dstat(0)).mean(), sd1 = (dg.dstat(0)).sd();
  double m2 = (dg.dstat(1)).mean(), sd2 = (dg.dstat(1)).sd();
  double xmin = 0.0, xmax = 0.0;
  double xrange = set_xrange(xmin, xmax, opt, dg);
  double dx = xrange / (double)RDdef;
  fout << "%理論的感度曲線 群1:コード=";
  if ((opt.g(0)).empty()) fout << opt.g(1) << " 以外\n";
  else fout << opt.g(0) << "\n";
  double x = xmin;
  for (uint i = 0; i <= RDdef; ++i, x += dx)
  {
    double sn = dg.fsn((x - m1) / sd1);
    fout << x << " " << sn << "\n";
  }
  fout << NGPdelm2 << "\n";
  fout << "%理論的特異度曲線 群2:コード=";
  if ((opt.g(1)).empty()) fout << opt.g(0) << " 以外\n";
  else fout << opt.g(1) << "\n";
  x = xmin;
  for (uint i = 0; i <= RDdef; ++i, x += dx)
  {
    double sp = dg.fsp((x - m2) / sd2);
    fout << x << " " << sp << "\n";
  }
}

int dpplotclass::get(const char* c_fname, option& opt, idd& id, ldata& datalst)
// データを読み込み累積する
{
  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("実行時オプションエラー！");
  nc = (int)(id.nc());          // ケース数変更用
  ldata_itr ditr = datalst.begin();
  while (ditr != datalst.end()) // リストの要素を削除するのでforではなくwhile
  {
    int idx = opt.search_g(ditr->sdata(0));
    double fdata = ditr->fdata(1);
    if (idx < 0 || fdata == BLANKf) // 欠損データはリストから削除する
    { ditr = datalst.erase(ditr);  --nc; }
    else
    {
      ditr->assign(0, (double)idx); // 群コードを群インデックスに変える
      ++ditr;
    }
  }
  if (nc < (int)(id.nc())) id.nc(nc); // ケース数を減らす
  return nc;
}

void dpplotclass::put(option& opt, idd& id, ldata& datalst) const
// 集計結果を出力する
{
  constexpr int LENline = 74;   // 区切り線の長さ
  datalst.sort(data_less(1));   // 項目データ順にソートする
  diag dg(opt.pi());            // 診断率オブジェクト
  dg.cal_dpa(opt, datalst);
  if (dg.resdiag() == reserr)
  {
    if (opt.res("res"))
    {
      put_resid(cout, id);
      if (opt.res("stat"))
      {
        cout << "\n";
        putline(cout, LENline);
        put_stat(cout, opt, dg);
        putline(cout, LENline);
      }
      cout << "※DP-pot分析不可能！\n\n" << endl;
    }
    return;
  }
  if (opt.res("res")) put_result(cout, opt, id, datalst, dg);
  if (opt.ngpname() != "") plot_ngp(opt, id, datalst, dg);
}

int main(int argc, char* argv[]) // メインルーチン
{
  dpplotclass dpplot;           // メインクラス変数
  option opt(argv[0]);          // オプション
  idd id;                       // データ識別情報
  ldata datalst;                // データリスト
  try
  {
    opt.setdef();               // オプションの初期設定
    bool fok = true;            // ファイル処理フラッグ
    int nc = 0;                 // 読み込んだケース数
    for (int i = 1; i <= argc; ++i)
    {
      if (i >= argc)
      {
        if (fok) dpplot.get(nullptr, opt, id, datalst); // 標準入力から読み込む
      }
      else if (*(argv[i]) == CHRopt)
      { opt.set(argv[i] + 1);  fok = true; }
      else
      { dpplot.get(argv[i], opt, id, datalst);  fok = false; }
    }
    dpplot.put(opt, id, datalst);
  }
  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;
}
