// cra.C Version 7.44  Copyleft (c) 杉本典夫
// Time-stamp: <2025-10-30 17:42:30 tac>
//
// 曲線回帰分析を行う。
//
// 使用ライブラリ:libdans.so, libfun.so
// コンパイル法 :g++ -O -s -W cra cra.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 =            // プログラム使用法表示用
{
  "cra  Version 7.44            Copyleft(Vo\\oV) 2025 杉本典夫",
  "書式:cra  [<オプション> …]  [<データファイル> …]",
  "機能:最初の項目を説明変数、次の項目を目的変数として曲線回帰分析をするよ。",
  "<オプション>",
  "  -cc[=<信頼係数>]:信頼係数を指定するよ。(0≦信頼係数<1)",
  "    -cc:-cc=0.95(オプション省略時)",
  "  -f=<ファイル名>:指定されたファイルからオプションを読み込むよ。",
  "  -h[=<コマンド名>]:オンラインマニュアルを表示するよ。(-h:-h=man)",
  "  -n=[<開始次数>][-[<終了次数>]][,…]:説明変数の次数を指定するよ。",
  "  -ngp={c|cl|data|pl}[,…]:Ngraph用出力データを指定するよ。",
  "    c:回帰曲線  cl:信頼限界  data:個々のデータ  pl:予測限界",
  "  -ngpname[=<出力ファイル名>]:Ngraph用データの出力ファイル名を指定するよ。",
  "    オプション省略時:Ngraph用データ出力せず  -ngpname:-ngpname=cra_ngp.dat",
  "  -res=[-][{def|all]|reg|stat}][,…]:解析結果の出力項目を指定するよ。",
  "    -:-を付けた項目を出力しない(-res=-:全ての項目を出力しない)",
  "    def:オプション省略時と同じ項目を出力する  all:全ての項目を出力する",
  "    reg:回帰分析結果  stat:基礎統計量",
  "  -xoffset[=<オフセット値>]:説明変数のオフセット値を指定するよ。",
  "  -xrange=<最小値>-<最小値>:X軸の最小値と最大値を指定するよ。"
};

const string Title = "=== 回帰曲線(regression curve) ==="; // 表題
const string NGPNAMEdef = "cra_ngp.dat"; // 省略時Ngraphファイル名
constexpr int RDdef = 50;       // グラフの座標範囲の分割数省略時値

enum ngptype                    // Ngraph用出力データの種類
{
  ngpnon,                       // 出力なし
  ngpc,                         // 回帰曲線
  ngpcl,                        // 信頼限界
  ngpdata,                      // 個々の生データ
  ngppl,                        // 予測限界
  ngpetc                        // その他
};

using listint = list<int>;      // 整数リストの短名定義
using listint_itr = listint::iterator; // 反復子の短名定義
using listint_citr = listint::const_iterator; // const反復子の短名定義

class option                    // オプション集合体
{
  string v_path;                // プログラムのパス
  string v_pname;               // プログラム名
  double v_cc;                  // 信頼係数
  list<pairint> v_nplst;        // 説明変数の次数ペアリスト
  listint v_nlst;               // 説明変数の次数リスト
  string v_ngpname;             // Ngraphファイル名
  lstr ngplst;                  // Ngraph出力リスト
  mapsb ngpmap;                 // Ngraph出力の有無マップ
  lstr reslst;                  // 分析結果出力リスト
  mapsb resmap;                 // 分析結果出力の有無マップ
  bool v_xoffset;               // 説明変数のオフセットフラッグ
  double v_x0;                  // 説明変数のオフセット値
  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 cc() const             // 信頼係数を返す
  { return v_cc; }
  list<pairint>& nplst()        // 説明変数の次数ペアリストを返す
  { return v_nplst; }
  listint& nlst()               // 説明変数の次数リストを返す
  { return v_nlst; }
  const string& ngpname() const // Ngraphデータファイル名を返す
  { return v_ngpname; }
  bool ngp(const string& code); // Ngraph出力フラッグを返す
  bool res(const string& code); // 分析結果出力フラッグを返す
  bool xoffset() const          // 説明変数のオフセットフラッグを返す
  { return v_xoffset; }
  double x0() const             // 説明変数のオフセット値を返す
  { return v_x0; }
  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 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_nplst(), v_nlst(), v_ngpname(), ngplst(), ngpmap(), reslst(), resmap(), v_xoffset(false), v_x0(0.0), xrange(BLANKf, 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); // プログラム名削除
}

bool option::ngp(const string& code) // Ngraph出力フラッグを返す
{
  mapsb_citr ngpitr = ngpmap.find(code);
  if (ngpitr == ngpmap.end()) return false; // code無し！
  return ngpmap[code];                      // true/falseを返す
}

bool option::res(const string& code) // 分析結果出力フラッグを返す
{
  mapsb_citr resitr = resmap.find(code);
  if (resitr == 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] == "f")       // オプション指定ファイル
  { err = fgetopt(*this, opt[1], true); }
  else if (opt[0] == "h")       // オンラインヘルプ表示
  { err = help(pname(), USAGE, opt[1]); }
  else if (opt[0] == "n")       // 説明変数の次数ペア指定
  {
    if (getrangelst(v_nplst, opt[1], CHRdelm) < 1) 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] == "res")     // 分析結果の出力指定
  {
    if (strsplit(reslst, opt[1], CHRdelm) < 1) err = badopt(pname(), optstr, USAGE);
  }
  else if (opt[0] == "xoffset") // 説明変数のオフセット値指定
  {
    v_xoffset = true;
    v_x0 = opt[1].empty() ? 0.0 : strtodf(opt[1]);
  }
  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オプションのセット
  if (v_nplst.size() < 1)       // 次数ペア指定なし
  {
    v_nlst.push_back(1);         // 1次式にする
    return 0;
  }
  bool nerr = false;
  list<pairint>::const_iterator ipitr = v_nplst.begin();
  while (ipitr != v_nplst.end() && nerr == false)
  {
    int nstart = ipitr->first, nend = ipitr->second;
    if (nstart < 1 || nend < 1)
    { nerr = true; }
    else
    {
      int nn = abs(nend - nstart) + 1;
      int inc = (nstart < nend) ? 1 : -1;
      int ni = nstart;
      for (int i = 0; i < nn; ++i)
      {
        v_nlst.push_back(ni);
        ni += inc;
      }
      ++ipitr;
    }
  }
  if (nerr) return disperr(pname(), "0以下の次数が指定されてるよ！");
  for (listint::const_iterator iitr = v_nlst.begin(); iitr != v_nlst.end(); ++iitr)
  {
    int ni = *iitr;
    listint::const_iterator jitr = iitr;
    ++jitr;
    while (jitr != v_nlst.end() && nerr == false)
    {
      if (ni == *jitr) nerr = true;
      else ++jitr;
    }
  }
  if (nerr) return disperr(pname(), "次数が重複してるよ！");
  return 0;                     // エラーなし
}

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

int option::set_res()           // -resオプションのセット
{
  resmap["res"] = resmap["def"] = resmap["reg"] = resmap["stat"] = true;
  if (reslst.size() > 0)
  {
    resmap["def"] = resmap["reg"] = resmap["stat"] = false;
    for (const string& resstr : reslst)
    {
      mapsb_citr resitr = resmap.find(resstr);
      if (resstr == "-")        // 全て出力せず
      { resmap["res"] = false; }
      else if (resstr[0] == '-' && resstr.size() > 1) // 特定の項目を出力せず
      {
        resitr = resmap.find(resstr.substr(1));
        if (resitr == resmap.end())
        { return disperr(pname(), "-res=" + resstr + "オプション指定が不適当！"); }
        else
        { resmap[resstr.substr(1)] = false; }
      }
      else if (resstr == "def" || resstr == "all")
      { resmap["def"]= resmap["reg"] = resmap["stat"] = true; }
      else if (resitr == resmap.end())
      { return disperr(pname(), "-res=" + resstr + "オプション指定が不適当！"); }
      else
      { resmap[resstr] = true; }
    }
  }
  return 0;                     // エラーなし
}

class varclass                 // 変数クラス
{
private:
  pairdbl xrange;               // Xの最小値と最大値
  statcor<double> statc;        // 相関基礎統計量
  mat<double> smat;            // 積和行列
  int v_nx;                     // 説明変数の数
  int v_np;                     // パラメータ数
  int v_dfr;                    // 残差自由度
  int v_dft;                    // 全体自由度
  double v_sb;                  // 回帰平方和
  double v_syy;                 // 全体平方和
  double v_vr;                  // 残差分散
  double v_aic;                 // AIC
  restype v_resvar;             // 処理結果の種類
public:
  varclass(option& opt);        // コンストラクタ
  ~varclass()                   // デストラクタ
  {}
  double xmin() const           // Xの最小値を返す
  { return xrange.first; }
  double xmax() const           // Xの最大値を返す
  { return xrange.second; }
  int nx() const                // 説明変数の数を返す
  { return v_nx; }
  int np() const                // パラメータ数を返す
  { return v_np; }
  int n() const                 // 例数を返す
  { return statc.n(); }
  double mx() const             // Xの平均値を返す
  { return statc.mx(); }
  double vxx() const            // Xの分散を返す
  { return statc.vxx(); }
  double my() const             // Yの平均値を返す
  { return statc.my(); }
  double vyy() const            // Yの分散を返す
  { return statc.vyy(); }
  double b(int idx) const       // 偏回帰係数を返す
  { return (idx <= np()) ? smat(idx, np()) : 0.0; }
  double vb(int idx) const      // 平方和を返す
  { return (idx <= np()) ? vr() * smat(idx, idx) : 0.0; }
  int dfb() const               // 回帰自由度を返す
  { return v_nx; }
  int dfr() const               // 残差自由度を返す
  { return v_dfr; }
  int dft() const               // 全体自由度を返す
  { return v_dft; }
  double sb() const             // 回帰平方和を返す
  { return v_sb; }
  double sr() const             // 残差平方和を返す
  { return smat(np(), np()); }
  double syy() const            // 全体平方和を返す
  { return v_syy; }
  double vr() const             // 残差分散を返す
  { return v_vr; }
  double r2() const             // 重寄与率を返す
  { return (syy() > 0.0) ? sb() / syy() : 0.0;; }
  double r2a() const            // 自由度調整済み重寄与率を返す
  { return (syy() > 0.0) ? fabs(1.0 - vr() * (double)dft() / syy()) : 0.0; }
  double aic() const            // AICを返す
  { return v_aic; }
  restype resvar() const        // 処理結果の種類を返す
  { return v_resvar; }
  void cal(option& opt, ldata& datalst); // データを集計して各種統計量を求める
  double calestimate(option& opt, double x) const;
  // Yの推定値を計算して返す
  double calcr(option& opt, double x, double t, ngptype ltype) const;
  // Yの限界幅を計算して返す
};

varclass::varclass(option& opt) : xrange(FDmax, -FDmax), statc(), smat(), v_nx(0), v_np(0), v_dfr(0), v_dft(0), v_sb(0.0), v_syy(0.0), v_vr(0.0), v_aic(0.0), v_resvar(resnon)
// コンストラクタ
{
  v_nx = (int)((opt.nlst()).size()); // 説明変数の数
  size_t dim = ((opt.xoffset()) ? 1 : 2) + v_nx;
  v_np = (int)dim - 1;
  smat.resize(dim, dim);
}

void varclass::cal(option& opt, ldata& datalst)
// データを集計して各種統計量を求める
{
  size_t dim = smat.r();
  vecdbl fdatavec(dim, 0.0);
  int offset = 1;
  if (opt.xoffset()) offset = 0;
  else fdatavec[0] = 1.0;
  ldata_itr ditr = datalst.begin();
  while (ditr != datalst.end())
  {
    double x = ditr->fdata(0), y = ditr->fdata(1);
    if (statc.sigma(x, y) < 1) // 欠損データはリストから削除する
    { ditr = datalst.erase(ditr); }
    else
    {
      if (x < xrange.first) xrange.first = x;
      if (x > xrange.second) xrange.second = x;
      if (opt.xoffset()) x -= opt.x0();
      listint_citr iitr = (opt.nlst()).begin();
      for (int i = 0; i < nx(); ++i, ++iitr) fdatavec[offset + i] = pow(x, *iitr);
      fdatavec[np()] = y;
      for (size_t i = 0; i < dim; ++i)
      {
        x = fdatavec[i];
        for (size_t j = i; j < dim; ++j) smat.add(i, j, x * fdatavec[j]);
      }
      ++ditr;
    }
  }
  for (int i = 0; i < np(); ++i)
  {
    for (size_t j = i + 1; j < dim; ++j) smat.assign(j, i, smat(i, j));
  }
  if (statc.n() < np())         // 回帰計算不可能
  {
    v_resvar = reserr;
    return;
  }
  for (int i = 0; i < np(); ++i) smat.sweep(i);
  if (smat(np(), np()) < FDeps) smat.assign(np(), np(), 0.0);
  if (opt.xoffset())
  {
    offset = 1;  v_dft = statc.n();
    v_syy = statc.syy() + (double)(statc.n()) * statc.my() * statc.my();
  }
  else
  {
    offset = 0;  v_dft = statc.n() - 1;
    v_syy = statc.syy();
  }
  v_dfr = v_dft - nx();
  v_vr = (v_dfr > 0) ? smat(np(), np()) / (double)v_dfr : 0.0;
  v_sb = v_syy - smat(np(), np());
  if (v_sb < FDeps) v_sb = 0.0;
  v_aic = (smat(np(), np()) > 0.0) ? (double)(statc.n()) * log(smat(np(), np()) / (double)(statc.n())) + (double)(2 * (nx() + 2 - offset)) : 0.0;
  v_resvar = resok;
}

double varclass::calestimate(option& opt, double x) const
// Yの推定値を計算して返す
{
  if (opt.xoffset()) x -= opt.x0();
  double yest = b(0);
  listint_citr iitr = (opt.nlst()).begin();
  if (opt.xoffset())
  {
    yest = b(0) * pow(x, *iitr);
    ++iitr;
  }
  for (int i = 1; i < np(); ++i, ++iitr) yest += b(i) * pow(x, *iitr);
  return yest;                  // 推定値を返す
}

double varclass::calcr(option& opt, double x, double t, ngptype ltype) const
// Yの限界幅を計算して返す
{
  vec<double> xv(smat.r());
  int offset = 1;
  if (opt.xoffset()) { x -= opt.x0();  offset = 0; }
  else xv.assign(0, 1.0);
  listint_citr iitr = (opt.nlst()).begin();
  for (int i = offset; i < np(); ++i, ++iitr) xv.assign(i, pow(x, (double)(*iitr)));
  xv.assign(np(), 0.0);
  vec<double> tv(xv.trn() * smat);
  double w = inp(tv, xv);
  if (ltype == ngppl) w += 1.0; // 予測限界幅
  return t * sqrt(v_vr * w);    // 限界幅を返す
}

class craclass                  // craメインクラス
{
private:
  string fname;                 // ファイル名
  idd id;                       // データ識別情報
  ldata datalst;                // データリスト
  void put_resid(ostream& fout);
  // 解析結果のタイトルとデータ識別情報を出力する
  void put_reg(ostream& fout, option& opt, const varclass& var);
  // 曲線回帰分析結果を出力する
  void plot_ngraph(option& opt, const varclass& var);
  // Ngraph用データ出力
  void plot_ngpid(ofstream& fout, option& opt);
  // Ngraph用データのタイトルとデータ識別情報を出力する
  void plot_ngpcurve(ofstream& fout, option& opt, const varclass& var);
  // Ngraph用に回帰曲線を出力する
  void plot_ngpdata(ofstream& fout);
  // Ngraph用に個々のデータを出力する
  void plot_ngplim(ofstream& fout, option& opt, ngptype ltype, const varclass& var);
  // Ngraph用に信頼・予測限界を出力する
public:
  craclass() : fname(), id(), datalst() // コンストラクタ
  {}
  ~craclass()                   // デストラクタ
  {}
  int get(const char* c_fname, option& opt); // データを読み込み累積する
  void put(option& opt);                     // 集計結果を出力する
};

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

void craclass::put_reg(ostream& fout, option& opt, const varclass& var)
// 曲線回帰分析結果を出力する
{
  fout << "\n・曲線回帰式";
  if (opt.xoffset()) fout << "  切片無:x=(x-" << opt.x0() << ")";
  fout << "\n";
  constexpr int Llen = 74;      // 区切り線の長さ
  fout << "変数    偏回帰係数      標準誤差";
  fout << setw(5) << opt.cc() * 100.0 << "%絶対精度  信頼区間下限  信頼区間上限\n";
  putline(fout, Llen);
  listint_citr iitr = (opt.nlst()).begin();
  for (int i = 0; i < var.np(); ++i)
  {
    double bi = var.b(i), vbi = var.vb(i);
    double ci = calmcci(mcfis, opt.cc(), 1, 1, var.dfr(), 1, vbi);
    if (i == 0)
    {
      if (opt.xoffset())
      {
        fout << "x^" << left << setw(2) << *iitr;
        ++iitr;
      }
      else
      { fout << "切片"; }
    }
    else
    {
      fout << "x^" << left << setw(2) << *iitr;
      ++iitr;
    }
    fout << right << setw(14) << bi << setw(14) << sqrt(vbi);
    fout << setw(14) << ci << setw(14) << bi - ci;
    fout << setw(14) << bi + ci << "\n";
  }
  putline(fout, Llen);
  fout << "重寄与率  :R2乗=" << var.r2();
  fout << "  自由度調整済重寄与率  :R'2乗=" << var.r2a() << "\n";
  fout << "重相関係数:R   =" << sqrt(var.r2()) << "  自由度調整済重相関係数:R'   =" << sqrt(var.r2a()) << "\n";
  fout << "AIC(赤池の情報量基準)=" << var.aic() << "\n\n";
  putatab(fout, "回帰", var.sb(), var.dfb(), var.sr(), var.dfr(), var.syy(), var.dft());
}

void craclass::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("data"))          // 個々のデータ
  {
    fout << "-s1 -dmark,00\n";
    plot_ngpid(fout, opt);
    plot_ngpdata(fout);
    delml = true;
  }
  if (opt.ngp("c"))             // 回帰曲線
  {
    if (delml)
    { fout << NGPdelm << "\n"; }
    else
    {
      fout << "-s1 -dline\n";
      delml = true;
    }
    plot_ngpid(fout, opt);
    plot_ngpcurve(fout, opt, var);
  }
  if (opt.ngp("cl"))            // 回帰曲線の信頼限界
  {
    if (delml)
    { fout << NGPdelm << "\n"; }
    else
    {
      fout << "-s1 -dline\n";
      delml = true;
    }
    plot_ngpid(fout, opt);
    plot_ngplim(fout, opt, ngpcl, var);
  }
  if (opt.ngp("pl"))            // 回帰曲線の予測限界
  {
    if (delml) fout << NGPdelm << "\n";
    else fout << "-s1 -dline\n";
    plot_ngpid(fout, opt);
    plot_ngplim(fout, opt, ngppl, var);
  }
}

void craclass::plot_ngpid(ofstream& fout, option& opt)
// Ngraph用データのタイトルとデータ識別情報を出力する
{
  fout << "%";  puttitle(fout, Title);
  fout << "%\n";
  id.put_title(fout, formngp);
  id.put_stra(fout, formngp);
  fout << "%説明変数x:" << id.item(0) << "\n";
  fout << "%目的変数y:" << id.item(1) << "\n";
  fout << "%説明変数の次数:";
  for (listint_citr iitr = (opt.nlst()).begin(); iitr != (opt.nlst()).end(); ++iitr) fout << *iitr << " ";
  if (opt.xoffset()) fout << " 切片無:x=(x-" << opt.x0() << ")";
  fout << "\n";
}

void craclass::plot_ngpdata(ofstream& fout)
// Ngraph用に個々のデータを出力する
{
  fout << "%生データ\n";
  for (auto& cdatai : datalst)
  { fout << cdatai.fdata(0) << " " << cdatai.fdata(1) << "\n"; }
}

void craclass::plot_ngpcurve(ofstream& fout, option& opt, const varclass& var)
// Ngraph用に回帰曲線を出力する
{
  fout << "%回帰曲線\n";
  double xmin = var.xmin(), xmax = var.xmax(), xmg = 0.5;
  if (opt.xmin() != BLANKf && opt.xmax() != BLANKf)
  { xmin = opt.xmin();  xmax = opt.xmax();  xmg = 0.0; }
  double xrange = xmax - xmin;
  double xstart = xmin - xrange * xmg, xend = xmax + xrange * xmg;
  double dx = (xend - xstart) / (double)RDdef;
  double x = xstart, y = 0.0;
  for (int i = 0; i <= RDdef; ++i, x += dx)
  {
    y = var.calestimate(opt, x);
    fout << x << " " << y << "\n";
  }
}

void craclass::plot_ngplim(ofstream& fout, option& opt, ngptype ltype, const varclass& var)
// Ngraph用に信頼・予測限界を出力する
{
  fout << "%回帰曲線の" << opt.cc() * 100;
  if (ltype == ngpcl) fout << "%信頼限界上限\n";
  else fout << "%予測限界上限\n";
  double t = ptot(1.0 - opt.cc(), var.dfr());
  double xmin = var.xmin(), xmax = var.xmax(), xmg = 0.5;
  if (opt.xmin() != BLANKf && opt.xmax() != BLANKf)
  { xmin = opt.xmin();  xmax = opt.xmax();  xmg = 0.0; }
  double xrange = xmax - xmin;
  double xstart = xmin - xrange * xmg, xend = xmax + xrange * xmg;
  double dx = (xend - xstart) / (double)RDdef;
  double x = xstart, y = 0.0;
  for (int i = 0; i <= RDdef; ++i, x += dx)
  {
    y = var.calestimate(opt, x) + var.calcr(opt, x, t, ltype);
    fout << x << " " << y << "\n";
  }
  fout << NGPdelm2 << "\n";
  fout << "%回帰曲線の" << opt.cc() * 100;
  if (ltype == ngpcl) fout << "%信頼限界下限\n";
  else fout << "%予測限界下限\n";
  x = xstart;
  for (int i = 0; i <= RDdef; ++i, x += dx)
  {
    y = var.calestimate(opt, x) - var.calcr(opt, x, t, ltype);
    fout << x << " " << y << "\n";
  }
}

int craclass::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 craclass::put(option& opt) // 集計結果を出力する
{
  constexpr int LENline = 69;   // 区切り線の長さ
  varclass var(opt);            // 変数オブジェクト
  var.cal(opt, datalst);        // 各種の計算実行
  if (opt.res("res"))           // 計算結果出力
  {
    put_resid(cout);
    if (var.n() < 1)
    { cout << "※有効なデータなし！\n"; }
    else if (opt.res("stat"))
    {
      cout << "\n・変数の基礎統計量\n";
      putline(cout, LENline);
      cout << "x:";  putmsd(cout, var.n(), var.mx(), var.vxx());
      cout << "y:";  putmsd(cout, var.n(), var.my(), var.vyy());
      putline(cout, LENline);
    }
    if (var.resvar() == reserr) cout << "※回帰曲線計算不可能！\n";
    else if (opt.res("reg")) put_reg(cout, opt, var);
    cout << "\n" << endl;
  }
  if (var.resvar() != reserr && opt.ngpname() != "") plot_ngraph(opt, var);
  // Ngraph用データ出力
}

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