// dpcv.C Version 7.43  Copyleft (c) 杉本典夫
// Time-stamp: <2025-10-29 09:13:27 tac>
//
// 最初の項目を群項目、次の項目を集計項目としてDP-plot分析の交差検証を行う。
//
// 使用ライブラリ:libdans.so, libfun.so
// コンパイル法  :g++ -O -s -W dpcv docv.C -I/usr/local/include
//               -L/usr/local/lib -ldans -lfun

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

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

const vecstr USAGE =            // プログラム使用法表示用
{
  "dpcv  Version 7.43           Copyleft(Vo\\oV) 2025 杉本典夫",
  "書式:dpcv [<オプション> …] [<データファイル> …]",
  "機能:最初の項目を群項目、次の項目を集計項目として",
  "     DP-plot分析の交差検証を行います。",
  "<オプション>",
  "  -f=<ファイル>:指定したファイルからオプションを読み込みます。",
  "  -g{1|2}[=<群コード>]:群1または群2のコードを指定します。",
  "  -h[=<コマンド名>]:オンラインマニュアルを表示します。(-h:-h=man)",
  "  -pi[=<群1の事前確率>]:群1の事前確率を指定します。(0<事前確率<1)",
  "  -r=<データ数>:検証群のデータ数を指定します。",
  "  -roc[={ac|cm|rp|yp|<境界値>}]:ROC分析用境界値の指定。(-roc:-roc=rp)"
};

constexpr int Nlim = 5;         // 境界値数
constexpr int Nh = Nlim + 1;    // パラメータ数
constexpr double PIdef = 0.5;   // 省略時事前確率

using vecstatv = vector<statv<double>>; // 基礎統計量ベクトルの別名定義
using lditr = list<ldata_itr>;  // データリスト反復子リストの別名定義
using pairvec = pair<vecdbl, vecint>; // ベクトルペアの別名定義
using ldpair = list<pairvec>;   // ベクトルペアリストの別名定義
using ldpair_itr = ldpair::iterator; // 反復子の別名定義
using ldpair_citr = ldpair::const_iterator; // const反復子の別名定義

double fsn(double z, bool uporder) // 理論的感度を返す
{
  return (uporder) ? 1.0 - (double)(nedp(z)) : (double)(nedp(z));
}

double fsp(double z, bool uporder) // 理論的特異度を返す
{
  return (uporder) ? nedp(z) : 1.0 - nedp(z);
}

double fac(double sn, double sp, double pi1) // 理論的正診率を返す関数
{
  return pi1 * sn + (1.0 - pi1) * sp;
}

double fppv(double sn, double sp, double pi1) // 理論的PPVを返す関数
{
  double ppv = (sn < FDeps || pi1 < FDeps) ? 1.0 : 1.0 / (1.0 + (1.0 - pi1) / pi1 * (1.0 - sp) / sn);
  if (ppv < FDeps) ppv = 0.0;
  return (ppv > 1.0) ? 1.0 : ppv;
}

double fnpv(double sn, double sp, double pi1) // 理論的NPVを返す関数
{
  double npv = (sp < FDeps || (1.0 - pi1) < FDeps) ? 1.0 : 1.0 / (1.0 + pi1 / (1.0 - pi1) * (1.0 - sn) / sp);
  if (npv < FDeps) npv = 0.0;
  return (npv > 1.0) ? 1.0 : npv;
}

enum roctype                    // ROC分析用境界値の種類
{
  rocnon,                       // ROC分析結果非出力
  rocac,                        // 正診率最大点
  roccm,                        // 最大χ2乗法
  rocrp,                        // ROC曲線最適点(R-point)
  rocyp,                        // Youden index(感度+特異度-1)最大点(Y-point)
  roclim,                       // 境界値(上下限は自動判定)
  rocetc                        // その他
};

class option                    // オプションクラス
{
private:
  string v_path;                // プログラムのパス
  string v_pname;               // プログラム名
  int v_r;                      // 検証群のデータ数
  vecstr gvec;                  // 群コード
  double v_pi;                  // 群1の事前確率
  roctype v_roc;                // ROC分析用境界値の種類
  double v_lim;                 // 境界値
public:
  option(const char* c_pname = nullptr); // コンストラクタ
  ~option()                              // デストラクタ
  {}
  const string& path() const    // プログラムのパスを返す
  { return v_path; }
  const string& pname() const   // プログラム名を返す
  { return v_pname; }
  int r() const                 // 検証群データ数を返す
  { return v_r; }
  const string& g(int idx) const // 群コードを返す
  { return gvec[idx]; }
  double pi() const             // 群1の事前確率を返す
  { return v_pi; }
  roctype roc() const           // ROC分析結果用境界値の種類を返す
  { return v_roc; }
  double lim() const            // 境界値を返す
  { return v_lim; }
  void setdef()                 // 省略時オプション指定を設定する
  { getoptdef(*this); }
  void set(const string& optstr, char devc = OPTdev);
  // オプション指定によってオプションを設定する
  int check(size_t ni);         // オプションのチェックとセット
};

option::option(const char* c_pname) : v_path(), v_pname(), v_r(1), gvec(2), v_pi(BLANKf), v_roc(rocrp), v_lim(BLANKf)
// コンストラクタ
{
  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); // プログラム名削除
  gvec[0] = "1";
}

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] == "f")            // オプション指定ファイル
  { err = fgetopt(*this, opt[1], true); }
  else if (opt[0] == "g1")      // 群1のコード指定
  { gvec[0] = opt[1]; }
  else if (opt[0] == "g2")      // 群2のコード指定
  { gvec[1] = opt[1]; }
  else if (opt[0] == "h")       // オンラインヘルプ表示
  { err = help(pname(), USAGE, opt[1]); }
  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] == "r")       // 検証群データ数の指定
  {
    char* ptr = nullptr;
    if (opt[1].empty()) v_r = 1;
    else v_r = (int)(strtol(opt[1].c_str(), &ptr, 10));
    if (v_r < 1) err = badopt(pname(), optstr, USAGE, 1.0, (double)INT_MAX);
  }
  else if (opt[0] == "roc")     // ROC分析結果の出力と境界値の種類の指定
  {
    if (opt[1].empty() || opt[1] == "rp") v_roc = rocrp;
    else if (opt[1] == "ac") v_roc = rocac;
    else if (opt[1] == "cm") v_roc = roccm;
    else if (opt[1] == "yp") v_roc = rocyp;
    else
    {
      v_roc = roclim;
      v_lim = strtodf(opt[1]);
      if (v_lim == BLANKf) err = badopt(pname(), optstr, USAGE);
    }
  }
  else                          // 間違ったオプション指定
  { err = nonopt(pname(), optstr, USAGE); }
  if (err > 0) throw help_except("ヘルプ表示");
  if (err < 0) throw option_error(string("オプション指定違反→") + optstr);
}

int option::check(size_t ni)    // オプションのセットとチェック
{
  if (ni < 2) return disperr(pname(), "項目が足りません！");
  if (gvec[0].empty() && gvec[1].empty()) return disperr(pname(), "群コードが指定されていません！");
  return 0;                     // エラーなし
}

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 varclass                  // 変数クラス
{
private:
  int v_ncv;                    // 検証回数
  int v_nlim;                   // 境界値数
  int v_nh;                     // パラメータ数
  vecint ndpvec;                // 検証結果例数ベクトル
  vecdbl hvec;                  // パラメータベクトル
  double v_pi10;                // 群1の事前確率(全例)
  double v_pi1;                 // 群1の事前確率
  bool v_uporder;               // 昇順フラッグ(m1>=m2の時true)
  ldpair v_cvdpairlst;          // 交差検証法の結果リスト
  int cal_dp(option& opt, lditr& lditrlstt);
  // DP-plot分析を行って総例数を返す
  double cal_acpoint(double lim0, double m1, double v1, double m2, double v2) const;
  // AC-pointを返す
  double cal_pnpoint(double lim0, double m1, double sd1, double m2, double sd2) const;
  // PN-pointを返す
  double cal_roc(double& auc, int n1, int n2, option& opt, lditr& lditrlstt) const;
  // ROC分析の境界を返す
  void cv_dp(lditr& lditrlstv);
  // DP-plot分析結果でデータリストを検証して正判別例数を合計する
  void set_cvdlst(int vno);     // 検証結果をリストに追加する
  int cal_dpcv(option& opt, ldata& datalst, int r, int k);
  // 交差検証を行って検証回数を返す
public:
  varclass(int nlim0, int nh0) : v_ncv(0), v_nlim(nlim0), v_nh(nh0), ndpvec(nlim0 * 4, 0), hvec(nh0, 0.0), v_pi1(0.0), v_uporder(true), v_cvdpairlst()
  // デフォールトコンストラクタ
  {}
  ~varclass()                   // デストラクタ
  {}
  int ncv() const               // 検証回数を返す
  { return v_ncv; }
  int nlim() const              // 境界値数を返す
  { return v_nlim; }
  int nh() const                // パラメータ数を返す
  { return v_nh; }
  double pi10() const           // 群1の事前確率を返す(全例)
  { return v_pi10; }
  double pi1() const            // 群1の事前確率を返す
  { return v_pi1; }
  bool uporder() const          // 昇順フラッグを返す
  { return v_uporder; }
  const ldpair& cvdpairlst() const // 交差検証法の結果リストを返す
  { return v_cvdpairlst; }
  void cal(option& opt, ldata& datalst); // データを集計して変数をセットする
};

int varclass::cal_dp(option& opt, lditr& lditrlstt)
// DP-plot分析を行って総例数を返す
{
  vecstatv statvec(2);          // 基礎統計量ベクトル
  for (ldata_itr ditr : lditrlstt) statvec[(size_t)(ditr->fdata(0))].sigma(ditr->fdata(1)); // 群ごとに統計量を集計する
  int n1 = statvec[0].n(), n2 = statvec[1].n();
  int nt = n1 + n2;
  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();
  double sdt = sd1 + sd2;
  double auc = 0.0;
  v_uporder = (m1 >= m2) ? true : false;
  v_pi1 = (opt.pi() == BLANKf) ? (double)n1 / (double)nt : opt.pi();
  double lim0 = (m1 + m2) * 0.5;  // 判別分析の境界値
  hvec[0] = (fabs(sdt) < FDeps) ? lim0 : (m1 * sd2 + m2 * sd1) / sdt;
  hvec[1] = cal_acpoint(lim0, m1, v1, m2, v2);
  hvec[2] = cal_pnpoint(lim0, m1, sd1, m2, sd2);
  hvec[3] = lim0;
  hvec[4] = cal_roc(auc, n1, n2, opt, lditrlstt);
  hvec[5] = auc;
  return n1 + n2;               // 総例数を返す
}

double varclass::cal_acpoint(double lim0, double m1, double v1, double m2, double v2) const
// AC-pointを返す
{
  if (sgndiff(m1, m2) == 0) return lim0; // 計算不可能、初期値を返す
  double vd = v2 - v1;
  if (sgndiff(v2, v1) == 0)     // 2群の分散が等しい時
  {
    return (m1 + m2) * 0.5 + v1 * log((1.0 - pi1()) / pi1()) / (m1 - m2);
  }
  if (v1 < FDeps || v2 < FDeps) return lim0; // 計算不可能、初期値を返す
  double md = m1 * v2 - m2 * v1;
  double sd1 = sqrt(v1), sd2 = sqrt(v2);
  double tmp = md * md - vd * (m1 * m1 * v2 - m2 * m2 * v1 + 2.0 * v1 * v2 * log((1.0 - v_pi1) / v_pi1 * sd1 / sd2));
  double limm = (md - sqrt(tmp)) / vd, limp = (md + sqrt(tmp)) / vd;
  return (uporder()) ? limm : limp; // 理論値を返す
}

double varclass::cal_pnpoint(double lim0, double m1, double sd1, double m2, double sd2) const
// PN-pointを返す
{
  if (sd1 < FDeps || sd2 < FDeps) return lim0; // 初期値を返す
  constexpr double Rerror = 1.0e-6;            // 収束限界
  constexpr double Delta = Rerror * 10.0; // 刻み幅
  constexpr double SDx = 5.0;   // 範囲を決めるためのSDの倍率
  double lim = lim0;
  double z1 = (lim - m1) / sd1, z2 = (lim - m2) / sd2;
  double sn = fsn(z1, uporder()), sp = fsp(z2, uporder());
  double pr = fppv(sn, sp, pi1()) - fnpv(sn, sp, pi1());
  while (fabs(pr) >= Rerror)
  {
    double dlim = Delta * lim;
    double dlim1 = dlim / sd1, dlim2 = dlim / sd2;
    double snd = fsn(z1 + dlim1, uporder()), spd = fsp(z2 + dlim2, uporder());
    double prd = fppv(snd, spd, pi1()) - fnpv(snd, spd, pi1());
    double g = (prd - pr) / dlim;
    if (fabs(g) < FDeps)
    { pr = 0.0; }
    else
    {
      lim -= pr / g;
      z1 = (lim - m1) / sd1;  z2 = (lim - m2) / sd2;
      sn = fsn(z1, uporder());  sp = fsp(z2, uporder());
      pr = fppv(sn, sp, pi1()) - fnpv(sn, sp, pi1());
    }
  }
  double lmin = m2 - sd1 * SDx, lmax = m1 + sd2 * SDx;
  if (uporder() == false) { lmin = m1 - sd1 * SDx;  lmax = m2 + sd2 * SDx; }
  if (lim < lmin || lmax < lim) return lim0; // 範囲外の時は初期値を返す
  return lim;                   // 範囲内
}

double varclass::cal_roc(double& auc, int n1, int n2, option& opt, lditr& lditrlstt) const
// ROC分析の境界を返す
{
  int nt = n1 + n2, tp = n1, dfn = 0, tn = 0, dtn = 0;
  int tp_max = 0, tn_max = n2;
  if (uporder() == false) lditrlstt.reverse(); // 群1の平均値の方が小さい時
  roctype roc = opt.roc();
  if (roc == rocnon) roc = rocrp; // ROC分析結果非出力時も分析だけはしておく
  double ac_max = 0.0, csqr_max = 0.0, yi_max = 0.0, dist_min = FDmax;
  double w = 0.0, lim = BLANKf;
  double fd = 0.0, fdpre = BLANKf;
  for (ldata_itr ditr : lditrlstt)
  {
    fd = ditr->fdata(1);
    if (fd != fdpre)
    {
      tp -= dfn;  tn += dtn;
      double sn = (double)tp / (double)n1;
      double sp = (double)tn / (double)n2;
      double ac = (double)(tp + tn) / (double)nt;
      w += (double)(tp * dtn) + 0.5 * (double)(dfn * dtn);
      double csqr = 0.0, yi = 0.0, x = 0.0, y = 0.0, dist = 0.0;
      switch (roc)
      {
      case rocac :              // 正診率最大
        if (ac > ac_max)
        { ac_max = ac;  lim = fd; }
        break;
      case roccm : csqr = csqr22(tn, n2 - tn, n1 - tp, tp); // 最大χ2乗法
        if (csqr > csqr_max)
        { csqr_max = csqr;  lim = fd; }
        break;
      case rocyp : yi = sn + sp - 1.0; // Youden index(Y-point)
        if (yi > yi_max)
        { yi_max = yi;  lim = fd; }
        break;
      default : x = 1.0 - sp;  y = 1.0 - sn; // ROC曲線最適点(R-point)
        dist = x * x + y * y;
        if (dist < dist_min)
        { dist_min = dist;  lim = fd; }
        break;
      }
      fdpre = fd;  dfn = dtn = 0;
    }
    if ((int)(ditr->fdata(0)) == 0) ++dfn;
    else ++dtn;
  }
  if (roc == roclim) lim = opt.lim();
  else if (lim == BLANKf) lim = fd; // 境界値がサーチできていない時
  w += 0.5 * (double)(dfn * dtn);
  auc = (n1 > 0 && n2 > 0) ? w / (double)n1 / (double)n2 : 1.0;
  return lim;                   // 境界値を返す
}

void varclass::cv_dp(lditr& lditrlstv)
// DP-plot分析結果でデータリストを検証して正判別例数を合計する
{
  ndpvec.assign(ndpvec.size(), 0);
  for (ldata_itr ditr : lditrlstv)
  {
    int gidx = (int)(ditr->fdata(0));
    double f = ditr->fdata(1);
    for (int i = 0; i < nlim(); ++i)
    {
      int j = i * 4;
      if (gidx == 0)              // 群1
      {
        ndpvec[j] += 1;
        if (uporder())
        {
          if (f >= hvec[i]) ndpvec[j + 1] += 1; // 正判別
        }
        else
        {
          if (f <= hvec[i]) ndpvec[j + 1] += 1; // 正判別
        }
      }
      else if (gidx == 1)         // 群2
      {
        ndpvec[j + 2] += 1;
        if (uporder())
        {
          if (f < hvec[i]) ndpvec[j + 3] += 1; // 正判別
        }
        else
        {
          if (f > hvec[i]) ndpvec[j + 3] += 1; // 正判別
        }
      }
    }
  }
}

void varclass::set_cvdlst(int vno) // 検証結果をリストに追加する
{
  vecdbl cvdhvec(1 + nh());     // 検証結果パラメータベクトル
  vecint cvdndpvec(nlim() * 6); // 検証結果例数ベクトル
  cvdhvec[0] = (double)((vno < 1) ? nh() : vno);
  for (int i = 0; i < nh(); ++i) cvdhvec[1 + i] = hvec[i];
  for (int i = 0; i < nlim(); ++i)
  {
    int j = i * 4, k = i * 6;
    cvdndpvec[k] = ndpvec[j];
    cvdndpvec[k + 1] = ndpvec[j + 1];
    cvdndpvec[k + 2] = ndpvec[j + 2];
    cvdndpvec[k + 3] = ndpvec[j + 3];
    cvdndpvec[k + 4] = ndpvec[j] + ndpvec[j + 2];
    cvdndpvec[k + 5] = ndpvec[j + 1] + ndpvec[j + 3];
  }
  v_cvdpairlst.push_back(make_pair(cvdhvec, cvdndpvec));
}

int varclass::cal_dpcv(option& opt, ldata& datalst, int r, int k)
// 交差検証を行って検証回数を返す
{
  lditr lditrlstt, lditrlstv; // トレーニングデータとバリデーションデータ反復子
  for (int i = 0; i < k; ++i)
  {
    int sidx = r * i, eidx = r * (i + 1) - 1;
    int j = 0;
    for (ldata_itr dpos = datalst.begin(); dpos != datalst.end(); ++dpos, ++j)
    {
      if (sidx <= j && j <= eidx) lditrlstv.push_back(dpos);
      else lditrlstt.push_back(dpos);
    }
    if (cal_dp(opt, lditrlstt) > 0)
    {
      cv_dp(lditrlstv);
      set_cvdlst(i + 1);
    }
    lditrlstt.clear();  lditrlstv.clear();
  }
  return k;                     // 検証回数を返す
}

void varclass::cal(option& opt, ldata& datalst)
// データを集計して変数をセットする
{
  lditr lditrlst;               // データリスト反復子リスト
  for (ldata_itr ditr = datalst.begin(); ditr != datalst.end(); ++ditr) lditrlst.push_back(ditr);
  int nt = cal_dp(opt, lditrlst);     // 全例のDP-plot分析(nt:総例数)
  v_pi10 = v_pi1;                     // 群1の確率(全例)
  int r = opt.r();                    // 検証群のデータ数
  v_ncv = 0;
  if (nt < 1 || r > nt - 2) return;   // 交差検証実行せず
  cv_dp(lditrlst);                    // 全例の検証
  set_cvdlst(0);                      // 全例の検証結果をリストに追加する
  lditrlst.clear();                   // 反復子リストを削除する
  int k = nt / r;                     // 検証群数
  if (nt % r > 0) ++k;                // 検証群数を切り上げる
  v_ncv = cal_dpcv(opt, datalst, r, k); // 交差検証を行って検証回数を保存する
}

class dpcvclass                 // dpcvメインクラス
{
private:
  string fname;                 // ファイル名
  idd id;                       // データ識別情報
  ldata datalst;                // データリスト
  int search_g(const option& opt, const string& gcode) const;
  // 群コードを検索して群インデックスを返す
  void put_dpcv(ostream& fout, const varclass& var);
  // 判別分析の交差検証結果をデータファイル形式で出力する
public:
  dpcvclass() : fname(), id(), datalst() // コンストラクタ
  {}
  ~dpcvclass()                  // デストラクタ
  {}
  int get(const char* c_fname, option& opt); // データを読み込み累積する
  void put(option& opt);                     // 集計結果を出力する
};

int dpcvclass::search_g(const option& opt, const string& gcode) const
// 群コードを検索して群インデックスを返す
{
  if (gcode == BLANKs) return -1; // 欠損データはヒットさせない
  int idx = -1;
  if ((opt.g(0)).empty())       // 群1のコード指定がない場合
  { idx = (gcode == opt.g(1)) ? 1 : 0; }
  else if (gcode == opt.g(0))   // 群1のコード
  { idx = 0; }
  else if ((opt.g(1)).empty())  // 群2のコード指定がない場合
  { idx = 1; }
  else if (gcode == opt.g(1))   // 群2のコード
  { idx = 1; }
  return idx;                   // -1:ヒットせず,0:群1,1:群2
}

void dpcvclass::put_dpcv(ostream& fout, const varclass& var)
// 判別分析の交差検証結果をデータファイル形式で出力する
{
  id.put_title(fout, formtag);
  id.put_stra(fout, formtag);
  idd tid(id.title(), 1 + var.nh() + var.nlim() * 6);
  tid.append_stra(id.item(0), atype, avalnon);
  tid.append_stra(id.item(1), atype, avalnon);
  tid.put_stra(fout, formtag);
  tid.assign_item(0, "パラメータ数/検証群番号");
  tid.assign_item(1, "SS-point");
  tid.assign_item(2, "AC-point");
  tid.assign_item(3, "PN-point");
  tid.assign_item(4, "判別分析の境界値");
  tid.assign_item(5, "ROC曲線の境界値");
  tid.assign_item(6, "ROC曲線のAUC");
  for (int i = 0; i < var.nlim(); ++i)
  {
    int j = 2 + var.nlim() + i * 6;
    tid.assign_item(j, tid.item(1 + i) + ":群1の例数");
    tid.assign_item(j + 1, tid.item(1 + i) + ":TP例数");
    tid.assign_item(j + 2, tid.item(1 + i) + ":群2の例数");
    tid.assign_item(j + 3, tid.item(1 + i) + ":FN例数");
    tid.assign_item(j + 4, tid.item(1 + i) + ":2群の合計例数");
    tid.assign_item(j + 5, tid.item(1 + i) + ":AC例数");
  }
  ostringstream ostr;
  ostr << "AC-point(群1の事前確率=" << var.pi10() << ")";
  tid.assign_item(2, ostr.str());
  ostr.str("");
  ostr << "PN-point(群1の事前確率=" << var.pi10() << ")";
  tid.assign_item(3, ostr.str());
  tid.put_item(fout, formtag);
  int prec0 = fout.precision(); // 現在の不動小数点精度(有効数字)を保存する
  for (ldpair_citr cvitr = (var.cvdpairlst()).begin(); cvitr != (var.cvdpairlst()).end(); ++cvitr)
  {
    fout << cvitr->first[0];
    fout.precision(16);
    for (size_t i = 1; i < cvitr->first.size(); ++i) fout << CHRdelm << cvitr->first[i];
    fout.precision(prec0);
    for (size_t i = 0; i < cvitr->second.size(); ++i) fout << CHRdelm << cvitr->second[i];
    fout << endl;
  }
}

int dpcvclass::get(const char* c_fname, option& opt) // データを読み込む
{
  if (c_fname == nullptr) fname.erase();
  else fname = c_fname;
  int nc = fgetdata(fname, id, datalst);
  if (nc < 0) throw runtime_error("実行時エラー！");
  if (opt.check(id.ni()) != 0) throw option_error("実行時オプションエラー！");
  nc = (int)(id.nc());          // ケース数変更用
  ldata_itr ditr = datalst.begin();
  while (ditr != datalst.end())
  {
    int idx = search_g(opt, 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 dpcvclass::put(option& opt) // 集計結果を出力する
{
  datalst.sort(data_less(1)); // データリストを項目データ順にソートする
  varclass var(Nlim, Nh);     // 変数オブジェクト
  var.cal(opt, datalst);
  if (var.ncv() > 0) put_dpcv(cout, var);
}

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