// drag.C Version 7.44  Copyleft (c) 杉本典夫
// Time-stamp: <2025-10-29 09:47:55 tac>
//
// グループ・ロジスティックモデルによる微分回帰分析を行う。
//
// 使用ライブラリ:libdans.so, libfun.so
// コンパイル法 :g++ -O -s -W drag drag.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 =            // プログラム使用法表示用
{
  "drag  Version 7.44          Copyleft(Vo\\oV) 2025 杉本典夫",
  "書式:drag  [<オプション> …]  [<データファイル> …]",
  "機能:1番目の項目を目的変数y、2番目の項目をyの変化量dy、3番目の項目を",
  "説明変数xとしてグループ・ロジスティックモデルによる微分回帰分析を行う。",
  "<オプション>", 
  "  -cc[=<信頼係数>]:信頼係数を指定する。(0≦信頼係数<1)",
  "    -cc:-cc=0.95(オプション省略時)",
  "  -di[=<平均感染期間>]:平均感染期間(1/γ)を指定する。(-di:-di=5)",
  "  -epsy[=<yの微小値>]:グラフ用yの微小値を指定する。(-epsy:-epsy=0.001)",
  "  -f=<ファイル名>:指定されたファイルからオプションを読み込む。",
  "  -h[=<コマンド名>]:オンラインマニュアルを表示する。(-h:-h=man)",
  "  -ngp={log|regd|rt|td}[,…]:Ngraph用出力データを指定する。",
  "    log:xとyのデータとロジスティック曲線",
  "    regd:xとdyのデータと微分ロジスティック曲線  rt:k値とRt曲線",
  "    td:K値と倍加時間曲線  オプション省略時:全て指定",
  "  -ngpname[=<出力ファイル名>]:Ngraph用データを出力するファイル名を指定する。",
  "    オプション省略時:データ出力せず  -ngpname:-ngpname=drag_ngp.dat",
  "  -range=<最小値>-<最大値>:Ngraph用Xの範囲を指定する。",
  "  -res=[-][{def|all|regc|reglog|stat}][,…]",
  "      :解析結果の出力項目を指定する。",
  "    -:-を付けた項目を出力しない(-res=-:全ての項目を出力しない)",
  "    def:オプション省略時と同じ項目を出力する  all:全ての項目を出力する",
  "    regc:yとdyの2次回帰分析  reglog:ロジスティック回帰分析  stat:基礎統計量",
  "    オプション省略時:データに合わせて自動判定",
  "  -simplex[={-|<最大反復回数>}]:滑降Simplex法の実行指定。(0≦最大反復回数)",
  "    -:実行せず  -simplex:-simplex=200(オプション省略時)",
  "  -xg=<開始値2>[,<開始値3>,…]:各グループのxの開始値を指定する。",
  "  -xu={log|real}:xの単位を指定する。",
  "    log:xを対数変換して計算  real:xを実数のまま計算(オプション省略時)"
};

const string Title = "=== 微分回帰分析(differential regression analysis) ==="; // 表題
const string NGPNAMEdef = "drag_ngp.dat"; // 省略時Ngraphファイル名
constexpr int Np = 2;             // 回帰式の次数
constexpr double DIdef = 5.0;     // 平均感染期間省略時値
constexpr int RMAXdef = 200;      // 最大計算回数省略時値
constexpr double EPSdef = 1.0e-5; // 収束誤差省略時値
constexpr int RDdef = 50;         // グラフの描画範囲分割数省略時値
constexpr double EPSYdef = 0.001; // yの微小値省略時値

enum xutype                     // xの単位の種類
{
  xupnon,                       // 単位なし
  xulog,                        // 対数
  xureal,                       // 実数
  xuetc                         // その他
};

double llogdis(double x, double m, double c, double y0)
// 対数ロジスティック分布関数
{
  if (x <= 0.0) return 0.0;
  return m / (1.0 + ((m - y0) / y0) * exp(-c * log(x)));
}

double llogden(double x, double m, double c, double y0)
// 対数ロジスティック分布密度関数
{
  if (x <= 0.0) return 0.0;
  double y = llogdis(x, m, c, y0);
  return c * y * (m - y) / m / x;
}

class option                    // オプション集合体
{
  string v_path;                // プログラムのパス
  string v_pname;               // プログラム名
  double v_cc;                  // 信頼係数
  double v_di;                  // 平均感染期間
  double v_epsy;                // yの微小値
  string v_ngpname;             // Ngraphファイル名
  lstr ngplst;                  // Ngraph出力リスト
  mapsb ngpmap;                 // Ngraph出力の有無マップ
  pairdbl xrange;               // Ngraph用Xの最小値と最大値
  lstr reslst;                  // 分析結果出力リスト
  mapsb resmap;                 // 分析結果出力の有無マップ
  bool v_simplexb;              // Simplex法実行フラッグ
  int v_simplexd;               // Simplex法最大計算回数
  lstr xglst;                   // xの開始値リスト
  vecdbl xgvec;                 // xの開始値ベクトル
  xutype v_xu;                  // xの単位
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; }
  double di() const             // 平均感染期間を返す
  { return v_di; }
  double epsy() const           // yの微小値を返す
  { return v_epsy; }
  bool ngp(const string& code); // Ngraph出力フラッグを返す
  const string& ngpname() const // Ngraphデータファイル名を返す
  { return v_ngpname; }
  size_t nxg() const            // 指定された開始値数を返す
  { return xgvec.size(); }
  bool res(const string& code); // 分析結果出力フラッグを返す
  double rmin() const           // Ngraph用Xの最小値を返す
  { return xrange.first; }
  double rmax() const           // Ngraph用Xの最大値を返す
  { return xrange.second; }
  bool simplexb() const         // Simplex法実行フラッグを返す
  { return v_simplexb; }
  int simplexd() const          // Simplex法最大計算回数を返す
  { return v_simplexd; }
  double xg(size_t idx) const   // グループの開始値を返す
  { return (idx < nxg()) ? xgvec[idx] : FDmax; }
  xutype xu() const             // xの単位を返す
  { return v_xu; }
  void setdef()                 // 省略時オプション指定を設定する
  { getoptdef(*this); }
  void set(const string& optstr, char devc = OPTdev);
  // オプション指定によってオプションを設定する
  int check(size_t ni);         // オプションのチェックとセット
  int set_ngp();                // -ngpオプションのセット
  int set_res();                // -resオプションのセット
};

option::option(const char* c_pname) : v_path(), v_pname(), v_cc(CCdef), v_di(DIdef), v_epsy(EPSYdef), v_ngpname(), ngplst(), ngpmap(), xrange(BLANKf, BLANKf), reslst(), resmap(), v_simplexb(true), v_simplexd(RMAXdef), xglst(), xgvec(), v_xu(xureal)
// コンストラクタ
{
  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::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] == "di")      // 平均感染期間
  {
    if (opt[1].empty()) v_di = DIdef;
    else v_di = strtodf(opt[1]);
    if (v_di <= 0.0) err = badopt(pname(), optstr, USAGE, 0.0, FDmax);
  }
  else if (opt[0] == "epsy")    // yの微小値
  {
    if (opt[1].empty()) v_epsy = EPSYdef;
    else v_epsy = strtodf(opt[1]);
    if (v_epsy < FDeps) err = badopt(pname(), optstr, USAGE, FDeps, FDmax);
  }
  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] == "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] == "range")  // Ngraph用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 if (opt[0] == "res")     // 分析結果の出力指定
  {
    if (strsplit(reslst, opt[1], CHRdelm) < 1) err = badopt(pname(), optstr, USAGE);
  }
  else if (opt[0] == "simplex") // Simplex法の実行指定
  {
    if (opt[1].empty())
    { v_simplexb = true;  v_simplexd = RMAXdef; }
    else if (opt[1] == "-")
    { v_simplexb = false;  v_simplexd = 0; }
    else
    {
      v_simplexb = true;
      char* ptr = nullptr;
      v_simplexd = strtol(opt[1].c_str(), &ptr, 0);
      if (v_simplexd < 0) err = badopt(pname(), optstr, USAGE, 0.0, (double)INT_MAX);
    }
  }
  else if (opt[0] == "xg")      // グループの開始値の指定
  {
    if (strsplit(xglst, opt[1], CHRdelm) < 1) err = badopt(pname(), optstr, USAGE);
  }
  else if (opt[0] == "xu")      // xの単位指定
  {
    if (opt[1].empty() || opt[1] == "real") v_xu = xureal;
    else if (opt[1] == "log") v_xu = xulog;
    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(size_t ni)    // オプションのセットとチェック
{
  if (ni != 3) return disperr(pname(), "項目は3個必要！");
  if (set_ngp() < 0) return -2; // -ngpオプションのセット
  if (set_res() < 0) return -3; // -resオプションのセット
  if (xglst.size() < 1) return 0; // グループの開始値指定なし
  xgvec.resize(xglst.size());
  xgvec.assign(xglst.size(), 0.0); // グループの開始値を0にする
  size_t i = 0;
  double xgpre = -FDmax;
  for (lstr_citr pos = xglst.begin(); pos != xglst.end(); ++pos, ++i)
  {
    xgvec[i] = strtodf(*pos);
    if (xgvec[i] <= xgpre || xgvec[i] == BLANKf) return disperr(pname(), "グループの開始値指定が不適切！");
    xgpre = xgvec[i];
  }
  return 0;                     // エラーなし
}

int option::set_ngp()           // -ngpオプションのセット
{
  ngpmap["log"] = ngpmap["regd"] = ngpmap["rt"] = ngpmap["td"] = true;
  if (ngplst.size() > 0)
  {
    if (v_ngpname == "") return disperr(pname(), "-ngpnameオプションの指定が必要！");
    ngpmap["log"] = ngpmap["regd"] = ngpmap["rt"] = ngpmap["td"] = 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["regc"] = resmap["reglog"] = resmap["stat"] = true;
  if (reslst.size() > 0)
  {
    resmap["def"] = resmap["regc"] = resmap["reglog"] = 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" || resstr == "all")
      {
        resmap["def"] = resmap["regc"] = resmap["reglog"] = resmap["stat"] = true;
      }
      else if (respos == resmap.end())
      { return disperr(pname(), "-res=" + resstr + "オプション指定が不適切！"); }
      else
      { resmap[resstr] = true; }
    }
  }
  return 0;                     // エラーなし
}

class varclass                  // 変数クラス
{
private:
  int v_np;                     // 回帰式の次数
  int v_ng;                     // グループ数
  pairdbl xrange;               // xの最小値と最大値
  pairdbl yrange;               // yの最小値と最大値
  statv<double> statvdy;        // dyの基礎統計量
  statcor<double> statcxy;      // xとyの相関基礎統計量
  statcor<double> statcydy;     // yとdyの相関基礎統計量
  mat<double> smat;             // 回帰計算用積和行列(回帰式は定数項なし)
  vecdbl mvec;                  // Mベクトル
  vecdbl cvec;                  // cベクトル
  vecdbl y0vec;                 // y0ベクトル
  restype v_resvar;             // 処理結果のフラッグ
  int cal_estmc(double& m, double& c, option& opt, ldata& datalstt, double xg);
  // yとdyの2次回帰分析によってmとcを推測する
  double cal_esty0(double m, double c, option& opt, ldata& datalstt, double xg);
  // ロジットの1次回帰分析によってy0を推測して返す
  restype cal_simplex(int& nr, vec<double>& hv, double m, double c, option& opt, ldata& datalstt, double xg);
  // Simplex法による非線形最小2乗計算
  void init_simplex(mat<double>& simplexmat, const vec<double>& hv);
  // シンプレックスの初期設定
  restype check_simplex(mat<double>& simplexmat, vec<double>& hlv, vec<double>& hmv, vec<double>& hhv, double& fl, double& fs, double& fh, int& lidx, int& hidx, int nr, double m, double c, option& opt, ldata& datalstt, double xg);
  // fiの計算と終了基準のチェック
  double cal_dyest(double y) const; // yの推定値を返す
  double cal_dycr(double y, double x, double vr, bool plim, option& opt) const;
  // dyの限界幅を返す
public:
  varclass(int np0, int ng0) : v_np(np0), v_ng(ng0), xrange(FDmax, -FDmax), yrange(FDmax, -FDmax), statvdy(), statcxy(), statcydy(), smat(np0 + 1, np0 + 1), mvec(ng0), cvec(ng0), y0vec(ng0), v_resvar(resnon)
  // デフォールトコンストラクタ
  {}
  ~varclass()                   // デストラクタ
  {}
  void setrange(ldata& datalst); // xとyの範囲をセットする
  int calreg(double& m, double& c, double& y0, int& nr, double& y00, double& c0, double& sry, option& opt, ldata& datalstt, double xg, int idx);
  // 各種の回帰分析を行って例数を返す
  int calmci(double& mlow, double& mup, double m, int dfr, double vr, option& opt) const;
  // Mの信頼区間を計算する(Newton法)
  double fq(const vec<double>& hv, double m, option& opt, ldata& datalstt, double xg) const;
  // 誤差関数
  int calxrange(double& xmin, double& xmax, double& dx, option& opt) const;
  // xの最小値と最大値とdxを返す
  double fnlogcum(double x, option& opt) const;
  // 累積ロジスティック関数値を返す
  double fnlogdcum(double x, option& opt) const;
  // 累積微分ロジスティック関数値を返す
  double fnlogd2cum(double x, option& opt) const;
  // 累積2階微分ロジスティック関数値を返す
  int np() const                // 回帰式の次数を返す(2)
  { return v_np; }
  int ng() const                // グループ数を返す
  { return v_ng; }
  double xmin() const           // xの最小値を返す
  { return xrange.first; }
  double xmax() const           // xの最大値を返す
  { return xrange.second; }
  double ymin() const           // yの最小値を返す
  { return yrange.first; }
  double ymax() const           // yの最大値を返す
  { return yrange.second; }
  int xyn() const               // xとyの例数を返す
  { return statcxy.n(); }
  double xymx() const           // xの平均値を返す
  { return statcxy.mx(); }
  double xyvxx() const          // xの不偏分散を返す
  { return statcxy.vxx(); }
  double xymy() const           // yの平均値を返す
  { return statcxy.my(); }
  double xyvyy() const          // yの不偏分散を返す
  { return statcxy.vyy(); }
  int dyn() const               // dyの例数を返す
  { return statvdy.n(); }
  double dymean() const         // dyの平均値を返す
  { return statvdy.mean(); }
  double dyvr() const           // dyの不偏分散を返す
  { return statvdy.vr(); }
  int ydyn() const              // yとdyの例数を返す
  { return statcydy.n(); }
  double ydymy() const          // yの平均値を返す
  { return statcydy.mx(); }
  double ydyvyy() const         // yの不偏分散を返す
  { return statcydy.vxx(); }
  double ydysyy() const         // yの平方和を返す
  { return statcydy.sxx(); }
  double ydymdy() const         // dyの平均値を返す
  { return statcydy.my(); }
  double ydyvdydy() const       // dyの不偏分散を返す
  { return statcydy.vyy(); }
  double ydysdydy() const       // dyの平方和を返す
  { return statcydy.syy(); }
  double ydysr() const          // 回帰式の残差を返す
  { return smat(np(), np()); }
  double ydyb(int idx) const // yとdyの偏回帰係数を返す
  { return (idx < np()) ? smat(idx, np()) : 0.0; }
  double ydys(int iidx, int jidx) const // yとdyの積和・平方和を返す
  { return (iidx < np() && jidx < np()) ? smat(iidx, jidx) : 0.0; }
  restype resvar() const       // 処理結果のフラッグを返す
  { return v_resvar; }
};

void varclass::setrange(ldata& datalst) // xとyの範囲をセットする
{
  for (auto& cdatai : datalst)
  {
    double y = cdatai.fdata(0), dy = cdatai.fdata(1), x = cdatai.fdata(2);
    if (y != BLANKf && dy != BLANKf && x != BLANKf)
    {
      if (x < xrange.first) xrange.first = x;
      if (x > xrange.second) xrange.second = x;
      if (y < yrange.first) yrange.first = y;
      if (y > yrange.second) yrange.second = y;
    }
  }
}

int varclass::calreg(double& m, double& c, double& y0, int& nr, double& y00, double& c0, double& sry, option& opt, ldata& datalstt, double xg, int idx)
// 各種の回帰分析を行って例数を返す
{
  int n = cal_estmc(m, c, opt, datalstt, xg); // mとcの推測
  mvec[idx] = m;  cvec[idx] = c;
  c0 = c;                       // 2次回帰式による推測値
  if (n < 1)                    // データなし
  {
    v_resvar = reserr;
    return n;
  }
  if (sgn(c) < 0 || sgn(m) <= 0) // cまたはM不適切
  {
    v_resvar = reserr;
    return n;
  }
  y0 = y00 = cal_esty0(m, c, opt, datalstt, xg);
  double e0 = (m - y0) / y0;
  vec<double> hv(2);            // y0とcをSimplex法のパラメターにする
  hv.assign(0, e0);  hv.assign(1, c);
  if (opt.simplexb())           // Simplex法実行
  {
    if ((v_resvar = cal_simplex(nr, hv, m, c, opt, datalstt, xg)) != reserr)
    {                           // Simplex法成功
      e0 = hv(0);  c = hv(1);
      y0 = m / (1.0 + e0);
    }
  }
  sry = fq(hv, m, opt, datalstt, xg); // 残差平方和
  if (n < 2 || sry < FDeps) sry = 0.0;
  mvec[idx] = m;  cvec[idx] = c;  y0vec[idx] = y0;
  return n;                      // 2次回帰分析の例数を返す
}

int varclass::calmci(double& mlow, double& mup, double m, int dfr, double vr, option& opt) const
// Mの信頼区間を計算する(Newton法)
{
  double t = ptot(1.0 - opt.cc(), dfr);
  double dy0 = 0.0, y = m * 1.1, dey = 1.0, slope = 1.0;
  double vpre = cal_dyest(m * 0.5);
  const double Rerror = EPSdef * vpre;
  constexpr double Delta = EPSdef;
  constexpr int Rmax = RMAXdef;
  int nr = 0;
  while (fabs((vpre = cal_dyest(y) + cal_dycr(y, t, vr, false, opt)) - dy0) >= Rerror && nr < Rmax)
  {
    dey = Delta * y;
    if (fabs(dey) < FDeps) break;
    slope = (cal_dyest(y + dey) + cal_dycr(y + dey, t, vr, false, opt) - vpre) / dey;
    if (fabs(slope) < FDeps) break;
    y += (dy0 - vpre) / slope;
    if (y <= 0.0) y = dey;
    ++nr;
  }
  mup = y;
  if (mup < m || FDmax < mup)   // 上限サーチできず
  {
    mup = mlow = BLANKf;
    return -1;                  // 信頼区間計算できず
  }
  y =  m - (mup - m);  nr = 0;
  while (fabs((vpre = cal_dyest(y) - cal_dycr(y, t, vr, false, opt)) - dy0) >= Rerror && nr < Rmax)
  {
    dey = Delta * y;
    if (fabs(dey) < FDeps) break;
    slope = (cal_dyest(y + dey) - cal_dycr(y + dey, t, vr, false, opt) - vpre) / dey;
    if (fabs(slope) < FDeps) break;
    y += (dy0 - vpre) / slope;
    if (y <= 0.0) y = dey;
    ++nr;
  }
  mlow = y;
  if (mlow < 0.0) mlow = 0.0;
  return nr;                    // 繰り返し数を返す
}

double varclass::fq(const vec<double>& hv, double m, option& opt, ldata& datalstt, double xg) const
// 誤差関数
{
  double q = 0.0;
  for (auto& cdatai : datalstt)
  {
    double y = cdatai.fdata(0), x = cdatai.fdata(2);
    if (y != BLANKf && x != BLANKf)
    {
      if (xg == BLANKf || x < xg) // グループデータのみ集計
      {
        if (opt.xu() == xulog) x = (x > 0.0) ? log(x) : 0.0;
        double ye = (sgn(hv(1)) == 0) ? m / (1.0 + hv(0)) : m / (1.0 + hv(0) * exp(-hv(1) * x));
        double yr = y - ye;     // hv(1)=c=0の時exp(-hv(1) * x)が1にならないため
        q += yr * yr;
      }
    }
  }
  return q;                     // 平方誤差合計を返す
}

int varclass::calxrange(double& xmin, double& xmax, double& dx, option& opt) const
// xの最小値と最大値とdxを返す
{
  int ng = opt.nxg() + 1;
  double m = mvec[0], c = cvec[0], y0 = y0vec[0];
  double ly0 = log((m - y0) / y0);
  double y = m * opt.epsy();    // y最小値
  xmin = (ly0 - log((m - y) / y)) / c;
  if (opt.xu() == xulog) xmin = exp(xmin);
  if (xmin > xrange.first) xmin = xrange.first;
  if (xmin < 0.0) xmin = 0.0;
  y = m * (1.0 - opt.epsy());   // y最大値
  xmax = (ly0 - log((m - y) / y)) / c;
  if (opt.xu() == xulog) xmax = exp(xmax);
  for (int i = 1; i < ng; ++i)
  {
    m = mvec[i];  c = cvec[i];  y0 = y0vec[i];
    ly0 = log((m - y0) / y0);
    y = m * (1.0 - opt.epsy());
    double x = (ly0 - log((m - y) / y)) / c;
    if (opt.xu() == xulog) x = exp(x);
    double xmaxi = opt.xg(i - 1) + x;
    if (xmax < xmaxi) xmax = xmaxi;
  }
  if (xmax < xrange.second) xmax = xrange.second;
  if (opt.rmin() != BLANKf) xmin = opt.rmin();
  if (opt.rmax() != BLANKf) xmax = opt.rmax();
  int rd = (int)(xmax - xmin);
  dx = 1.0;
  if (rd < RDdef)               // 分割数が省略時値未満の時
  {
    int r = (int)((double)RDdef / (double)rd + 0.9);
    rd *= r;  dx /= (double)r;
  }
  else if (rd > RDdef * 10)     // 分割数が省略時値の10倍以上の時
  {
    int r = (int)((double)rd / (double)(RDdef * 10) + 0.9);
    rd /= r;  dx *= (double)r;
  }
  return rd;                    // 分割数を返す
}

double varclass::fnlogcum(double x, option& opt) const
// 累積ロジスティック関数値を返す
{
  double y = 0.0, ye = 0.0, xg = 0.0;
  for (int i = 0; i < ng(); ++i)
  {
    double m = mvec[i], c = cvec[i], y0 = y0vec[i];
    if (m > 0.0)
    {
      double xt = x - xg;
      if (opt.xu() == xulog) ye = (xt > 0.0) ? llogdis(xt, m, c, y0) : 0.0;
      else ye = logdis(xt, m, c, y0);
      y += ye;
    }
    xg = opt.xg(i);
  }
  return y;
}

double varclass::fnlogdcum(double x, option& opt) const
// 累積微分ロジスティック関数値を返す
{
  double dy = 0.0, dye = 0.0, xg = 0.0;
  for (int i = 0; i < ng(); ++i)
  {
    double m = mvec[i], c = cvec[i], y0 = y0vec[i];
    if (m > 0.0)
    {
      double xt = x - xg;
      if (opt.xu() == xulog) dye = (xt > 0.0) ? llogden(xt, m, c, y0) : 0.0;
      else dye = logden(xt, m, c, y0);
      dy += dye;
    }
    xg = opt.xg(i);
  }
  return dy;
}

double varclass::fnlogd2cum(double x, option& opt) const
// 累積2階微分ロジスティック関数値を返す
{
  double d2y = 0.0, d2ye = 0.0, xg = 0.0;
  for (int i = 0; i < ng(); ++i)
  {
    double m = mvec[i], c = cvec[i], y0 = y0vec[i];
    double xt = x - xg;
    if (m > 0.0)
    {
      if (opt.xu() == xulog) d2ye = (xt > 0.0) ? (c * (m - 2.0 * llogdis(xt, m, c, y0)) / m - 1.0) / xt * llogden(xt, m, c, y0) : 0.0;
      else d2ye = c / m * (m - 2.0 * logdis(xt, m, c, y0)) * logden(xt, m, c, y0);
      d2y += d2ye;
    }
    xg = opt.xg(i);
  }
  return d2y;
}

int varclass::cal_estmc(double& m, double& c, option& opt, ldata& datalstt, double xg)
// yとdyの2次回帰分析によってmとcを推測する
{
  int d = np() + 1;
  size_t dyidx = np();
  m = yrange.second;  c = 0.0;
  smat = 0.0;  statvdy.clear();  statcxy.clear();  statcydy.clear();
  vecdbl fdatavec(d);
  for (auto& cdatai : datalstt)
  {
    double y = cdatai.fdata(0), dy = cdatai.fdata(1), x = cdatai.fdata(2);
    if (y != BLANKf && x != BLANKf)
    {
      statcxy.sigma(x, y);
      if (dy != BLANKf)
      {
        if (xg == BLANKf || x < xg) // グループデータのみ集計
        {
          statvdy.sigma(dy);              // dy累積
          if (opt.xu() == xulog) dy *= x; // 対数の時はx・dyにする
          statcydy.sigma(y, dy);
          double fd = y;
          for (int i = 0; i < np(); ++i)
          {
            fdatavec[i] = fd;
            fd *= y;
          }
          fdatavec[dyidx] = dy;
          for (int i = 0; i < d; ++i)
          {
            fd = fdatavec[i];
            for (int j = i; j < d; ++j) smat.add(i, j, fd * fdatavec[j]);
          }
        }
      }
    }
  }
  int n = statcydy.n();
  if (n < 1) return n;          // 2次回帰分析用データなし
  for (size_t i = 0; i < dyidx; ++i)
  {
    for (size_t j = i + 1; j <= dyidx; ++j) smat.assign(j, i, smat(i, j));
  }
  if (n < np()) return -1;      // yとdyの2次回帰分析不可能
  for (size_t i = 0; i < dyidx; ++i) smat.sweep(i);
  double b1 = smat(0, dyidx), b2 = smat(1, dyidx);
  if (n == np() || smat(dyidx, dyidx) < FDeps) smat.assign(dyidx, dyidx, 0.0);
  c = b1;  m = -b1 / b2;
  return n;                     // 2次回帰分析の例数を返す
}

double varclass::cal_esty0(double m, double c, option& opt, ldata& datalstt, double xg)
// ロジットの1次回帰分析によってy0を推測して返す
{
  statcor<double> statcxly;     // xとロジットの相関基礎統計量
  for (auto& cdatai : datalstt)
  {
    double y = cdatai.fdata(0), x = cdatai.fdata(2);
    if (y != BLANKf && x != BLANKf)
    {
      if (xg == BLANKf || x < xg) // グループデータのみ集計
      {
        if (opt.xu() == xulog) x = (x > 0.0) ? log(x) : 0.0;
        double ly = 0.0;
        if (y <= 0.0) ly = log(opt.epsy() / (1.0 - opt.epsy()));
        else if (m <= y) ly = log((1.0 - opt.epsy()) / opt.epsy());
        else ly = log(y / (m - y));
        statcxly.sigma(x, ly);
      }
    }
  }
  double y0 = yrange.first;
  if (statcxly.n() < 1) return y0; // ロジットがない時
  double ly0 = statcxly.my() - c * statcxly.mx();
  y0 = m / (1.0 + exp(-ly0));
  if (y0 >= m) return yrange.first; // 不適切なy0
  return y0;                    // y0推測値を返す
}

restype varclass::cal_simplex(int& nr, vec<double>& hv, double m, double c, option& opt, ldata& datalstt, double xg)
// Simplex法による非線形最小2乗計算
{
  int p = hv.dim();             // パラメーター数
  mat<double> simplexmat(p, p + 1); // シンプレックスを行列で表す
  vec<double> hlv(hv), hmv(hv), hhv(hv);
  double fl = 0.0, fs = 0.0, fh = 0.0; // 最小値、2番目に大きい値、最大値
  int lidx = 0, hidx = 0;
  restype ressimp = resnon;
  nr = 0;
  init_simplex(simplexmat, hv);
  while ((ressimp = check_simplex(simplexmat, hlv, hmv, hhv, fl, fs, fh, lidx, hidx, nr, m, c, opt, datalstt, xg)) == resnon)
  {
    ++nr;
    vec<double> hrv(hmv);
    hrv *= 2.0;
    hrv -= hhv;
    double fr = fq(hrv, m, opt, datalstt, xg);
    if (fr < fl)                // 折り返し成功
    {
      vec<double> hev(hrv);
      hev *= 2.0;
      hev -= hmv;
      double fe = fq(hev, m, opt, datalstt, xg);
      if (fe < fr)              // 拡大成功
      { for (int i = 0; i < p; ++i) simplexmat.assign(i, hidx, hev(i)); }
      else                      // 拡大失敗(行き過ぎ)
      { for (int i = 0; i < p; ++i) simplexmat.assign(i, hidx, hrv(i)); }
    }
    else if (fr <= fs)          // 斜め方向の折り返し
    { for (int i = 0; i < p; ++i) simplexmat.assign(i, hidx, hrv(i)); }
    else                        // 縮小失敗、Simplex移動
    {
      if (fr < fh) { fh = fr;  hhv = hrv; } // 折り返し後のSimplex中に最小点
      vec<double> hcv(hmv + hhv);
      hcv *= 0.5;
      double fc = fq(hcv, m, opt, datalstt, xg);
      if (fc < fh)              // 縮小成功
      { for (int i = 0; i < p; ++i) simplexmat.assign(i, hidx, hcv(i)); }
      else                      // 縮小失敗、Simplex移動
      {
        for (int j = 0; j <= p; ++j)
        {
          if (j != lidx)
          {
            for (int i = 0; i < p; ++i) simplexmat.assign(i, j, (simplexmat(i, j) + hlv(i)) * 0.5);
          }
        }
      }
    }
  }
  if (ressimp == reserr) return ressimp; // 発散
  hv = hlv;
  return ressimp;               // 収束または反復回数最大値
}

void varclass::init_simplex(mat<double>& simplexmat, const vec<double>& hv)
// シンプレックスの初期設定
{
  constexpr double Sqrt2 = sqrt(2.0); // コンパイル時定数
  constexpr double Rs = 0.1;          // シンプレックスの1辺の相対割合
  int p = hv.dim();
  double s = hv(0) * Rs;
  double v = (sqrt((double)(p + 1)) - 1.0 + (double)p) / (double)p / Sqrt2 * s;
  double w = (sqrt((double)(p + 1)) - 1.0) / (double)p / Sqrt2 * s;
  for (int i = 0; i < p; ++i) simplexmat.assign(i, 0, hv(i));
  for (int j = 1; j <= p; ++j)
  {
    for (int i = 0; i < p; ++i) simplexmat.assign(i, j, (i == j - 1) ? v + hv(i) : w + hv(i));
  }
}

restype varclass::check_simplex(mat<double>& simplexmat, vec<double>& hlv, vec<double>& hmv, vec<double>& hhv, double& fl, double& fs, double& fh, int& lidx, int& hidx, int nr, double m, double c, option& opt, ldata& datalstt, double xg)
// fiの計算と終了基準のチェック
{
  int p = hlv.dim();
  int eidx = p - 1;
  vecdbl fvec(p + 1);
  fl = FDmax;  fh = -FDmax;
  for (int j = 0; j <= p; ++j)
  {
    vec<double> hjv(simplexmat.sub(0, j, eidx, j));
    double fj = fq(hjv, m, opt, datalstt, xg);
    fvec[j] = fj;
    if (fj < fl) { fl = fj;  lidx = j; }
    if (fj > fh) { fh = fj;  hidx = j; }
  }
  if (fl > FDmax || isnan(fl)) return reserr; // 発散
  hlv = simplexmat.sub(0, lidx, eidx, lidx); // 最小値の時のhv
  hhv = simplexmat.sub(0, hidx, eidx, hidx); // 最大値の時のhv
  if (fabs(hlv(0)) > FDmax || isnan(hlv(0))) return reserr; // 発散
  double ferr = fabs(fh - fl);
  if (ferr < fabs(fl * EPSdef)) return resok; // 収束
  else if (opt.simplexd() > 0 && nr >= opt.simplexd()) return resrmax; // 反復回数最大値
  fs = fl;
  hmv = 0.0;
  for (int j = 0; j <= p; ++j)
  {
    if (j != hidx)
    {
      if (fvec[j] > fs) fs = fvec[j];
      hmv += simplexmat.sub(0, j, eidx, j);
    }
  }
  hmv *= 1.0 / (double)p;       // hhv以外のp個のhvの重心
  return resnon;                // 収束せず
}

double varclass::cal_dyest(double y) const // yの推定値を返す
{
  double yp = y;
  size_t dyidx = np();
  double vest = 0.0;
  for (size_t i = 0; i < dyidx; ++i)
  {
    vest += smat(i, dyidx) * yp;
    yp *= y;
  }
  return vest;                  // y推定値を返す
}

double varclass::cal_dycr(double y, double x, double vr, bool plim, option& opt) const
// dyの限界幅を返す
{
  if (opt.xu() == xulog && sgn(x) <= 0) return 0.0;
  size_t dyidx = np();
  vec<double> yv(dyidx + 1);
  double yp = y;
  for (size_t i = 0; i < dyidx; ++i)
  {
    yv.assign(i, yp);
    yp *= y;
  }
  yv.assign(dyidx, 0.0);
  vec<double> tv(yv.trn() * smat);
  double w = inp(tv, yv);
  if (plim) w += 1.0;           // 予測限界幅
  if (opt.xu() == xulog) x = log(x);
  return x * sqrt(vr * w);      // vの限界幅を返す
}

class dragclass                 // dragメインクラス
{
private:
  string fname;                 // ファイル名
  idd id;                       // データ識別情報
  ldata datalst;                // データリスト
  void put_resid(ostream& fout); // 分析結果のタイトルとデータ識別情報を出力する
  void put_stat(ostream& fout, option& opt, int ng, int gidx, double xgpre, double xg, const varclass& var) const;
  // 基礎統計量を出力する
  void put_regc(ostream& fout, option& opt, double& mt, double& mtlow, double& mtup, double m, double c, const varclass& var);
  // yとdyの2次回帰分析の結果を出力する
  void put_coef2(ostream& fout, option& opt, double& mt, double& mtlow, double& mtup, double m, double c, int dfr, double vr, const varclass& var);
  // 2次回帰式の偏回帰係数を出力する
  void put_reglog(ostream& fout, option& opt, double m, double c, double y0, int nr, double y00, double c0, double sry, ldata& datalstt, double xg, const varclass& var);
  // ロジスティック回帰分析結果を出力する
  void update_datalst(option& opt, ldata& datalstt, double m, double c, double y0, double xg);
  // 計算用データリストを更新する
  void plot_ngraph(option& opt, const varclass& var); // Ngraph用データ出力
  void plot_ngpid(ofstream& fout, option& opt);
  // Ngraph用データのタイトルとデータ識別情報を出力する
  void plot_ngpxy(ofstream& fout); // Ngraph用にxとyのデータを出力する
  void plot_ngplog(ofstream& fout, option& opt, const varclass& var) const;
  // Ngraph用にxとyのロジスティック曲線を出力する
  void plot_ngpxdy(ofstream& fout); // Ngraph用にxとdyのデータを出力する
  void plot_ngpregd(ofstream& fout, option& opt, const varclass& var) const;
  // Ngraph用にxとdyの微分ロジスティック曲線を出力する
  void plot_ngpkl(ofstream& fout, option& opt, const varclass& var) const;
  // Ngraph用にK値曲線を出力する
  void plot_ngptd(ofstream& fout, option& opt, const varclass& var) const;
  // Ngraph用に倍加時間曲線を出力する
  void plot_ngpks(ofstream& fout, option& opt, const varclass& var) const;
  // Ngraph用にk値曲線を出力する
  void plot_ngprt(ofstream& fout, option& opt, const varclass& var) const;
  // Ngraph用にRt曲線を出力する
public:
  dragclass() : fname(), id(), datalst() // コンストラクタ
  {}
  ~dragclass()                  // デストラクタ
  {}
  int get(const char* c_fname, option& opt); // データを読み込み累積する
  void put(option& opt);                     // 集計結果を出力する
};

void dragclass::put_resid(ostream& fout)
// 分析結果のタイトルとデータ識別情報を出力する
{
  puttitle(fout, Title);
  fout << "\n";
  id.put_title(fout);
  id.put_stra(fout);
  fout << "目的変数y  :" << id.item(0) << "\n";
  fout << "yの変化量dy:" << id.item(1) << "\n";
  fout << "説明変数x  :" << id.item(2) << endl;
}

void dragclass::put_stat(ostream& fout, option& opt, int ng, int gidx, double xgpre, double xg, const varclass& var) const
// 基礎統計量を出力する
{
  constexpr int LENline = 71;   // 区切り線の長さ
  fout << "\n○変数の基礎統計量";
  if (ng > 1) fout << ":グループ" << gidx + 1 << "(" << xgpre << "≦x＜" << xg << ")";
  fout << "\n";
  putline(fout, LENline);
  fout << "y  :";  putmsd(fout, var.xyn(), var.xymy(), var.xyvyy());
  fout << "dy :";  putmsd(fout, var.dyn(), var.dymean(), var.dyvr());
  if (opt.xu() == xulog)
  { fout << "xdy:";  putmsd(fout, var.ydyn(), var.ydymdy(), var.ydyvdydy()); }
  fout << "x  :";  putmsd(fout, var.xyn(), var.xymx(), var.xyvxx());
  putline(fout, LENline);
}

void dragclass::put_regc(ostream& fout, option& opt, double& mt, double& mtlow, double& mtup, double m, double c, const varclass& var)
// yとdyの2次回帰分析の結果を出力する
{
  fout << "\n・yとdyの2次回帰分析  モデル関数:";
  if (opt.xu() == xulog) fout << "x・dy=b1・y^1+b2・y^2\n";
  else fout << "dy=b1・y^1+b2・y^2\n";
  if (var.resvar() == reserr)
  {
    fout << "※計算エラー発生！" << endl;
    return;
  }
  int n = var.ydyn();
  size_t dyidx = var.np();
  double mdy = var.ydymdy();
  double st = var.ydysdydy() + (double)n * mdy * mdy;
  double sr = var.ydysr();
  int dft = n, dfb = 2;
  double sb = st - sr;
  int dfr = dft - dfb;
  double vr = (dfr > 0) ? sr / (double)dfr : 0.0;
  put_coef2(fout, opt, mt, mtlow, mtup, m, c, dfr, vr, var);
  double r2 = (st > 0.0) ? sb / st : 0.0;
  fout << "重寄与率(決定係数):R2乗=" << r2 << "  重相関係数:R=" << sqrt(r2) << "\n\n";
  putatab(fout, "回帰", sb, dfb, sr, dfr, st, dft);
}

void dragclass::put_coef2(ostream& fout, option& opt, double& mt, double& mtlow, double& mtup, double m, double c, int dfr, double vr, const varclass& var)
// 2次回帰式の偏回帰係数とパラメータを出力する
{
  fout << "記号    偏回帰係数      標準誤差";
  fout << setw(3) << opt.cc() * 100.0 << "%信頼区間：下限          上限\n";
  constexpr int Llen = 64;      // 区切り線の長さ
  putline(fout, Llen);
  size_t dyidx = var.np(), idx = 0, offset = 1;
  vecdbl blow(3, 0.0), bup(3, 0.0);
  for (size_t i = 0; i < dyidx; ++i)
  {
    double bi = var.ydyb(i), sii = var.ydys(i, i);
    double vbi = vr * sii;
    fout << "b" << left << setw(3) << offset + i;
    fout << right << setw(14) << bi << setw(14) << sqrt(vbi);
    double ci = calmcci(mcfis, opt.cc(), 1, 1, dfr, 1, vbi);
    blow[idx] = bi - ci;  bup[idx] = bi + ci;
    fout << setw(18) << blow[idx] << setw(14) << bup[idx] << "\n";
    ++idx;
  }
  putline(fout, Llen);
  double mlow = 0.0, mup = 0.0, clow = blow[0], cup = bup[0];
  var.calmci(mlow, mup, m, dfr, vr, opt);
  fout << "M   " << right << setw(14) << m;
  if (mup != BLANKf) fout << setw(32) << mlow << setw(14) << mup;
  fout << "\n";
  fout << "c   " << right << setw(14) << c;
  fout << setw(32) << clow << setw(14) << cup << "\n";
  if (opt.nxg() > 0)
  {
    mt += m;
    if (mup == BLANKf) { mtlow += m;  mtup += m; }
    else { mtlow += mlow;  mtup += mup; }
    fout << "Mt  " << right << setw(14) << mt;
    fout << setw(32) << mtlow << setw(14) << mtup << "\n";
  }
  putline(fout, Llen);
}

void dragclass::put_reglog(ostream& fout, option& opt, double m, double c, double y0, int nr, double y00, double c0, double sry, ldata& datalstt, double xg, const varclass& var)
// ロジスティック回帰分析結果を出力する
{
  fout << "\n・xとyのロジスティック回帰分析\n";
  fout << "  モデル関数:";
  if (opt.xu() == xulog) fout << "y=M/[1+{(M-y0)/y0}exp(-c・log(x))]=M/{1+exp(-ly0-c・log(x))}\n";
  else fout << "y=M/[1+{(M-y0)/y0}exp(-c・x)]=M/{1+exp(-ly0-c・x)}\n";
  fout << "  計算方法:";
  if (opt.simplexb())
  {
    fout << "非線形最小2乗法(滑降Simplex法)  計算回数=" << nr;
    if (var.resvar() == reserr)
    {
      fout << ":反復計算発散！" << endl;
      return;
    }
    else if (var.resvar() == resrmax)
    { fout << "(最大反復回数)\n"; }
    else
    { fout << "\n"; }
  }
  else
  { fout << "ロジットの最小2乗法\n"; }
  int n = var.xyn();
  double sr = sry;
  vec<double> hv(2);
  hv.assign(0, 1.0);  hv.assign(1, 0.0);
  double st = var.fq(hv, m, opt, datalstt, xg); // M=m,y0=m/2,c=0の時の残差をstにする
  double sb = st - sr;
  if (sb < FDeps) sb = 0.0;
  int dft = n - 1, dfb = 1;
  int dfr = dft - dfb;
  double ly00 = log(y00 / (m - y00));
  double ly0 = log(y0 / (m - y0));
  int llen = 74;
  putline(fout, llen);
  fout << "初期値:M=" << m << "(固定)  c=" << c0 << "  y0=" << y00 << "  ly0=" << ly00 << "\n";
  fout << "収束後:M=" << m << "(固定)  c=" << c << "  y0=" << y0 << "  ly0=" << ly0 << "\n";
  double y = m * 0.5;
  double ly = log(y / (m - y));
  double x = (ly - ly0) / c;
  double dy = logden(x, m, c, y0);
  if (opt.xu() == xulog) { dy = llogden(x, m, c, y0);  x = exp(x); }
  fout << "y=M/2=" << y << "の時のx=" << x << "  dy=" << dy << "\n";
  putline(fout, llen);
  double r2 = (sgn(st) > 0) ? sb / st : 0.0;
  fout << "重寄与率(決定係数):R2乗=" << r2 << "  重相関係数:R=" << sqrt(r2) << "\n\n";
  putatab(fout, "回帰", sb, dfb, sr, dfr, st, dft);
}

void dragclass::update_datalst(option& opt, ldata& datalstt, double m, double c, double y0, double xgt)
// 計算用データリストを更新する
{
  for (auto& cdatai : datalstt)
  {
    double y = cdatai.fdata(0), dy = cdatai.fdata(1), x = cdatai.fdata(2);
    if (x != BLANKf)
    {
      if (x < xgt)
      { cdatai.assign(2, BLANKf); }
      else
      {
        if (y != BLANKf)
        {
          double ye = (opt.xu() == xulog) ? llogdis(x, m, c, y0) : logdis(x, m, c, y0);
          cdatai.assign(0, y - ye);
        }
        if (dy != BLANKf)
        {
          double dye = (opt.xu() == xulog) ? llogden(x, m, c, y0) : logden(x, m, c, y0);
          cdatai.assign(1, dy - dye);
        }
        cdatai.assign(2, x - xgt);
      }
    }
  }
}

void dragclass::plot_ngraph(option& opt, const varclass& var)
// Ngraph用データ出力
{
  ofstream fout;
  if (ofopen(fout, (opt.ngpname()).c_str()) == false) return; // ファイルオープン失敗
  bool delml = false;
  if (opt.ngp("log"))
  {
    fout << "-s1 -dline\n";
    plot_ngpid(fout, opt);
    plot_ngpxy(fout);
    fout << NGPdelm << "\n";
    plot_ngplog(fout, opt, var);
    delml = true;
  }
  if (opt.ngp("regd"))
  {
    if (delml)
    { fout << NGPdelm << "\n"; }
    else
    {
      fout << "-s1 -dline\n";
      delml = true;
    }
    plot_ngpid(fout, opt);
    plot_ngpxdy(fout);
    fout << NGPdelm << "\n";
    plot_ngpregd(fout, opt, var);
  }
  if (opt.ngp("td"))
  {
    if (delml)
    { fout << NGPdelm << "\n"; }
    else
    {
      fout << "-s1 -dline\n";
      delml = true;
    }
    plot_ngpid(fout, opt);
    plot_ngpkl(fout, opt, var);
    fout << NGPdelm << "\n";
    plot_ngptd(fout, opt, var);
  }
  if (opt.ngp("rt"))
  {
    if (delml) fout << NGPdelm << "\n";
    else fout << "-s1 -dline\n";
    plot_ngpid(fout, opt);
    plot_ngpks(fout, opt, var);
    fout << NGPdelm << "\n";
    plot_ngprt(fout, opt, var);
  }
}

void dragclass::plot_ngpid(ofstream& fout, option& opt)
// Ngraph用データのタイトルとデータ識別情報を出力する
{
  fout << "%";  puttitle(fout, Title);
  fout << "%\n";
  id.put_title(fout, formngp);
  id.put_stra(fout, formngp);
  fout << "%目的変数y:" << id.item(0);
  fout << "\n%計算方法:";
  if (opt.simplexb()) fout << "非線形最小2乗法(滑降Simplex法)";
  fout << "\n%yの変化量dy:" << id.item(1) << "\n";
  fout << "%説明変数x:" << id.item(2) << endl;
}

void dragclass::plot_ngpxy(ofstream& fout)
// Ngraph用にxとyのデータを出力する
{
  fout << "%説明変数xと目的変数yのデータ\n";
  int dp = 1, i = 0;
  if (datalst.size() > RDdef * 10) // データ数が省略時分割数の10倍以上の時
  { dp = (int)((double)(datalst.size()) / (double)(RDdef * 10) + 0.9); }
  for (auto& cdatai : datalst)
  {
    double y = cdatai.fdata(0), x = cdatai.fdata(2);
    if (y != BLANKf && x != BLANKf)
    {
      if ((i % dp) == 0) { fout << x << " " << y << endl;  ++i; }
    }
  }
}

void dragclass::plot_ngplog(ofstream& fout, option& opt, const varclass& var) const
// Ngraph用にxとyのロジスティック曲線を出力する
{
  fout << "%説明変数xと目的変数yのロジスティック曲線\n";
  double xmin = 0.0, xmax = 1.0, dx = 1.0;
  int rd = var.calxrange(xmin, xmax, dx, opt);
  double x = xmin;
  for (int i = 0; i <= rd; ++i, x += dx)
  { fout << x << " " << var.fnlogcum(x, opt) << endl; }
}

void dragclass::plot_ngpxdy(ofstream& fout)
// Ngraph用にxとdyのデータを出力する
{
  fout << "%説明変数xと目的変数変化量dyのデータ\n";
  int dp = 1, i = 0;
  if (datalst.size() > RDdef * 10) // データ数が省略時分割数の10倍以上の時
  { dp = (int)((double)(datalst.size()) / (double)(RDdef * 10) + 0.9); }
  for (auto& cdatai : datalst)
  {
    double dy = cdatai.fdata(1), x = cdatai.fdata(2);
    if (dy != BLANKf && x != BLANKf)
    {
      if ((i % dp) == 0)
      {
        fout << x - 0.1 << " 0 " << x + 0.1 << " " << dy << endl;
        ++i;
      }
    }
  }
}

void dragclass::plot_ngpregd(ofstream& fout, option& opt, const varclass& var) const
// Ngraph用にxとdyの微分ロジスティック曲線を出力する
{
  fout << "%説明変数xと目的変数変化量dyの微分ロジスティック曲線\n";
  double xmin = 0.0, xmax = 1.0, dx = 1.0;
  int rd = var.calxrange(xmin, xmax, dx, opt);
  double x = xmin;
  for (int i = 0; i <= rd; ++i, x += dx)
  { fout << x << " " << var.fnlogdcum(x, opt) << endl; }
}

void dragclass::plot_ngpkl(ofstream& fout, option& opt, const varclass& var) const
// Ngraph用にK値曲線を出力する
{
  fout << "%目的変数yのK値曲線\n";
  double xmin = 0.0, xmax = 1.0, dx = 1.0;
  int rd = var.calxrange(xmin, xmax, dx, opt);
  double x = xmin;
  for (int i = 0; i <= rd; ++i, x += dx)
  {
    double kl = var.fnlogdcum(x, opt) / var.fnlogcum(x, opt);
    fout << x << " " << kl << endl;
  }
}

void dragclass::plot_ngptd(ofstream& fout, option& opt, const varclass& var) const
// Ngraph用に倍加時間曲線を出力する
{
  fout << "%目的変数yの倍加時間曲線\n";
  double xmin = 0.0, xmax = 1.0, dx = 1.0;
  int rd = var.calxrange(xmin, xmax, dx, opt);
  double x = xmin;
  for (int i = 0; i <= rd; ++i, x += dx)
  {
    double td = log(2.0) * var.fnlogcum(x, opt) / var.fnlogdcum(x, opt);
    fout << x << " " << td << endl;
  }
}

void dragclass::plot_ngpks(ofstream& fout, option& opt, const varclass& var) const
// Ngraph用にk値曲線を出力する
{
  fout << "%目的変数変化量dyのk値曲線\n";
  double xmin = 0.0, xmax = 1.0, dx = 1.0;
  int rd = var.calxrange(xmin, xmax, dx, opt);
  double x = xmin;
  for (int i = 0; i <= rd; ++i, x += dx)
  {
    double ks = var.fnlogd2cum(x, opt) / var.fnlogdcum(x, opt);
    fout << x << " " << ks << endl;
  }
}

void dragclass::plot_ngprt(ofstream& fout, option& opt, const varclass& var) const
// Ngraph用にRt曲線を出力する
{
  fout << "%目的変数変化量dyのRt曲線\n";
  double xmin = 0.0, xmax = 1.0, dx = 1.0;
  int rd = var.calxrange(xmin, xmax, dx, opt);
  double x = xmin;
  for (int i = 0; i <= rd; ++i, x += dx)
  {
    double rt = 1.0 + opt.di() * var.fnlogd2cum(x, opt) / var.fnlogdcum(x, opt);
    if (rt < 0.0) rt = 0.0;
    fout << x << " " << rt << endl;
  }
}

int dragclass::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("実行時オプションエラー！");
  return nc;
}

void dragclass::put(option& opt) // 集計結果を出力する
{
  if (opt.res("res")) put_resid(cout);
  int ng = opt.nxg() + 1;       // グループ数
  varclass var(Np, ng);         // 変数オブジェクト
  var.setrange(datalst);        // xとyの範囲をセットする
  ldata datalstt(datalst);      // 計算用データリスト
  int nt = 0;
  double mt = 0.0, mtlow = 0.0, mtup = 0.0;
  double xgpre = 0.0;
  for (int i = 0; i < ng; ++i)
  {
    int nr = 0;
    double y00 = 0.0, c0 = 0.0, m = 0.0, c = 0.0, y0 = 0.0, sry = 0.0;
    double xg = BLANKf, xgt = BLANKf;
    if (ng > 1)
    {
      xg = (i < ng - 1) ? opt.xg(i) : var.xmax() + 1.0;
      xgt = xg - xgpre;
    }
    int n = var.calreg(m, c, y0, nr, y00, c0, sry, opt, datalstt, xgt, i);
    nt += n;
    if (opt.res("res"))         // 計算結果出力
    {
      if (i > 0) cout << "\n";
      if (opt.res("stat")) put_stat(cout, opt, ng, i, xgpre, xg, var);
      if (n == 0)
      { cout << "※有効なデータなし！\n"; }
      else
      {
        if (opt.res("regc")) put_regc(cout, opt, mt, mtlow, mtup, m, c0, var);
        if (opt.res("reglog")) put_reglog(cout, opt, m, c, y0, nr, y00, c0, sry, datalstt, xgt, var);
      }
    }
    if (ng > 1)
    {
      update_datalst(opt, datalstt, m, c, y0, xgt);
      xgpre = xg;
    }
  }
  if (opt.res("res")) cout << "\n" << endl;
  if (nt > 0 && opt.ngpname() != "") plot_ngraph(opt, var);
}

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