// nlra.C Version 7.44  Copyleft (c) 杉本典夫
// Time-stamp: <2025-10-29 10:38:56 tac>
//
// 非線形最小二乗法による回帰分析
//
// 使用ライブラリ:libdans.so, libfun.so
// コンパイル法:g++ -O -s -W nlra nlra.C
//              -I/usr/local/include -L/usr/local/lib -ldans -lfun

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

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

const vecstr USAGE =            // プログラム使用法表示用
{
  "nlra  Version 7.44           Copyleft(Vo\\oV) 2025 杉本典夫",
  "書式:nlra  [<オプション> …]  [<データファイル> …]",
  "機能:最初の項目を説明変数、次の項目を目的変数として",
  "     非線形最小二乗法による回帰分析を行います。",
  "<オプション>",
  "  -eps=<収束誤差>:非線形最小二乗法の収束誤差を指定します。(0＜収束誤差)",
  "  -f=<ファイル名>:指定されたファイルからオプションを読み込みます。",
  "  -h[=<コマンド名>]:オンラインマニュアルを表示します。(-h:-h=man)",
  "  -h0=<初期値1>,…,<初期値p>:パラメーターの初期値を指定します。",
  "  -meth={newton|non|simplex|etc}:計算法を指定します。",
  "  -model={exp|log|log1|logc|etc}:モデルを指定します。",
  "  -ngp={data|func}[,…]:Ngraph用出力データを指定します。",
  "    data:説明変数と目的変数の生データ  func:モデル関数",
  "    オプション省略時:全て指定",
  "  -ngpname[=<出力ファイル名>]:Ngraph用データを出力するファイル名を指定します。",
  "    オプション省略時:データ出力せず  -ngpname:-ngpname=nlra_ngp.dat",
  "  -res=[-][{def|all|atab|matab|model|stat}][,…]:解析結果の出力項目指定。",
  "    -:-を付けた項目を出力しない(-res=-:全ての項目を出力しない)",
  "    def:オプション省略時と同じ項目を出力する  all:全ての項目を出力する",
  "    atab:分散分析表  matab:モデル固有の分散分析表  model:モデル関数",
  "    stat:基礎統計量  オプション省略時:データに合わせて自動判定",
  "  -rmax=<最大計算回数>:最大計算回数を指定します。(0≦最大計算回数)",
  "  -xrange=<最小値>-<最小値>:X軸の最小値と最大値を指定します。",
  "  -ym=<最小値>:yの最小値を指定します。(logcのみ)"
};

const string Title = "=== 非線形最小2乗法による曲線回帰分析 ===";
const string NGPNAMEdef = "nlra_ngp.dat"; // 省略時Ngraphファイル名
constexpr double EPSdef = 1.0e-4; // 収束誤差省略時値
constexpr int RMAXdef = 200;      // 最大計算回数省略時値
constexpr int RDdef = 50;         // グラフの描画範囲分割数省略時値
constexpr size_t Pdef = 2;        // 回帰式の次数

enum methtype                   // 非線形最小2乗法の種類
{
  methnon,                      // 指定なし
  methnewton,                   // Newton-Raphson法
  methsimplex,                  // Simplex法
  methetc                       // その他
};

enum modeltype                  // モデルの種類
{
  modelnon,                     // 指定なし
  modelexp,                     // 指数関数
  modellog,                     // ロジスティック関数
  modellog1,                    // ロジスティック関数(M=1固定)
  modellogc,                    // ロジスティック関数(最小値有)
  modeletc                      // その他
};

class option                    // オプション集合体
{
private:
  string v_path;                // プログラムのパス
  string v_pname;               // プログラム名
  double v_eps;                 // 非線形最小2乗法の収束誤差
  methtype v_meth;              // 非線形最小2乗法の種類
  modeltype v_model;            // モデルの種類
  int v_nh;                     // パラメータ数
  lstr h0lst;                   // 初期値リスト
  vecdbl h0vec;                 // 初期値ベクトル
  string v_ngpname;             // Ngraphファイル名
  lstr ngplst;                  // Ngraph出力リスト
  mapsb ngpmap;                 // Ngraph出力の有無マップ
  lstr reslst;                  // 分析結果出力リスト
  mapsb resmap;                 // 分析結果出力の有無マップ
  int v_rmax;                   // 最大計算回数
  pairdbl xrange;               // X軸の最小値と最大値
  double v_ym;                  // yの定数項
public:
  option(const char* c_pname = nullptr); // コンストラクタ
  ~option()                              // デストラクタ
  {}
  const string& path() const    // プログラムのパスを返す
  { return v_path; }
  const string& pname() const   // プログラム名を返す
  { return v_pname; }
  double eps() const            // 非線形最小2乗法の収束誤差を返す
  { return v_eps; }
  double h0(size_t idx) const   // 初期値を返す
  { return (idx < nh0()) ? h0vec[idx] : 0.0; }
  methtype meth() const         // 計算法の種類を返す
  { return v_meth; }
  const string methname() const; // 計算法の名称を返す
  modeltype model() const        // モデルの種類を返す
  { return v_model; }
  const string modelname() const; // モデルの名称を返す
  int nh() const                // パラメーター数を返す
  { return v_nh; }
  size_t nh0() const            // 指定された初期値数を返す
  { return h0vec.size(); }
  bool ngp(const string& code); // Ngraph出力フラッグを返す
  const string& ngpname() const // Ngraphデータファイル名を返す
  { return v_ngpname; }
  bool res(const string& code); // 分析結果出力フラッグを返す
  int rmax() const              // 最大計算回数を返す
  { return v_rmax; }
  double xmin() const           // X軸の最小値を返す
  { return xrange.first; }
  double xmax() const           // X軸の最大値を返す
  { return xrange.second; }
  double ym() const               // yの最小値を返す
  { return v_ym; }
  double ym(double ym0)         // yの最小値を設定する
  { return (v_ym = ym0); }
  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オプションのセット
  double fq(const vec<double>& hv, ldata& datalst) const;
  // 誤差関数
  double fqd(int hidx, const vec<double>& hv, ldata& datalst) const;
  // 誤差の偏導関数
  double fqd2(int hiidx, int hjidx, const vec<double>& hv, ldata& datalst) const;
  // 誤差の2階偏導関数
  double func(double x, const vec<double>& hv) const;
  // モデル関数
  double funcd(double x, int hidx, const vec<double>& hv) const;
  // モデル偏導関数
  double funcd2(double x, int hiidx, int hjidx, const vec<double>& hv) const;
  // モデル2階偏導関数
  double absmin(const vec<double>& hv) const;
  // 絶対値最小のパラメターの値を返す
};

option::option(const char* c_pname) : v_path(), v_pname(), v_eps(EPSdef), v_meth(methetc), v_model(modellog), v_nh(3), h0lst(), h0vec(), v_ngpname(), ngplst(), ngpmap(), reslst(), resmap(), v_rmax(RMAXdef), xrange(BLANKf, BLANKf), v_ym(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 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] == "eps")          // 非線形最小二乗法の収束誤差の指定
  {
    v_eps = opt[1].empty() ? EPSdef : strtodf(opt[1]);
    if (v_eps < 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] == "h0")      // 初期値の指定
  {
    if (strsplit(h0lst, opt[1], CHRdelm) < 1) err = badopt(pname(), optstr, USAGE);
  }
  else if (opt[0] == "meth")    // 計算法の指定
  {
    if (opt[1].empty() || opt[1] == "etc") v_meth = methetc;
    else if (opt[1] == "newton") v_meth = methnewton;
    else if (opt[1] == "non") v_meth = methnon;
    else if (opt[1] == "simplex") v_meth = methsimplex;
    else err = badopt(pname(), optstr, USAGE);
  }
  else if (opt[0] == "model")   // モデルの指定
  {
    if (opt[1].empty() || opt[1] == "log") { v_model = modellog;  v_nh = 3; }
    else if (opt[1] == "etc") { v_model = modeletc;  v_nh = 2; }
    else if (opt[1] == "exp") { v_model = modelexp;  v_nh = 2; }
    else if (opt[1] == "log1") { v_model = modellog1;  v_nh = 2; }
    else if (opt[1] == "logc") { v_model = modellogc;  v_nh = 3; }
    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] == "ym")      // yの最小値指定
  {
    v_ym = opt[1].empty() ? BLANKf : strtodf(opt[1]);
    if (v_ym == BLANKf) 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] == "rmax")    // 最大計算回数指定
  {
    char* ptr = nullptr;
    v_rmax = strtol(opt[1].c_str(), &ptr, 0);
    if (v_rmax < 0) err = badopt(pname(), optstr, USAGE, 0.0, (double)INT_MAX);
  }
  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 (h0lst.size() < 1) return 0; // 初期値指定なし
  if ((int)(h0lst.size()) > nh()) return disperr(pname(), "初期値の値が多すぎます！");
  h0vec.resize(nh(), 0.0);      // 初期値を0にする
  size_t i = 0;
  for (const string& h0str : h0lst)
  {
    h0vec[i] = strtodf(h0str);
    if (h0vec[i] == BLANKf) return disperr(pname(), "初期値指定が不適切です！");
    ++i;
  }
  return 0;                     // エラーなし
}

int option::set_ngp()           // -ngpオプションのセット
{
  ngpmap["data"] = ngpmap["func"] = true;
  if (ngplst.size() > 0)
  {
    if (v_ngpname == "") return disperr(pname(), "-ngpnameオプションの指定が必要！");
    ngpmap["data"] = ngpmap["func"] = 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["atab"] = resmap["model"] = resmap["stat"] = true;
  resmap["matab"] = false;
  if (reslst.size() > 0)
  {
    resmap["def"] = resmap["atab"] = resmap["model"] = 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") // 省略時出力
      {
        resmap["def"] = resmap["atab"] = resmap["model"] = resmap["stat"] = true;
      }
      else if (resstr == "all") // 全て出力
      {
        resmap["def"] = resmap["atab"] = resmap["matab"] = resmap["model"] = resmap["stat"] = true;
      }
      else if (respos == resmap.end())
      { return disperr(pname(), "-res=" + resstr + "オプション指定が不適切！"); }
      else
      { resmap[resstr] = true; }
    }
  }
  return 0;                     // エラーなし
}

const string option::methname() const // 計算法の名称を返す
{
  switch (meth())
  {
  case methnewton : return "Newton法";
  case methnon : return "なし";
  case methsimplex : return "滑降Simplex法";
  default : return "滑降Simplex+Newton法";
  }
}

const string option::modelname() const // モデルの名称を返す
{
  switch (model())
  {
  case modelexp : return "指数関数";
  case modellog : return "ロジスティック関数";
  case modellog1 : return "ロジスティック関数(M=1固定)";
  case modellogc : return "ロジスティック関数(最小値有)";
  case modelnon : return "なし";
  default : return "その他のモデル(1次関数)";
  }
}

double option::fq(const vec<double>& hv, ldata& datalst) const
// 誤差関数
{
  double q = 0.0;
  for (auto& cdatai : datalst)
  {
    double x = cdatai.fdata(0), y = cdatai.fdata(1);
    if (x != BLANKf && y != BLANKf)
    {
      double ei = y - func(x, hv);
      q += ei * ei;
    }
  }
  return q;                     // 平方誤差合計を返す
}

double option::fqd(int hidx, const vec<double>& hv, ldata& datalst) const
// 誤差の偏導関数
{
  double dq = 0.0;
  for (auto& cdatai : datalst)
  {
    double x = cdatai.fdata(0), y = cdatai.fdata(1);
    if (x != BLANKf && y != BLANKf)
    {
      double dei = y * funcd(x, hidx, hv) - func(x, hv) * funcd(x, hidx, hv);
      dq -= 2.0 * dei;
    }
  }
  return dq;                    // 誤差導関数値を返す
}

double option::fqd2(int hiidx, int hjidx, const vec<double>& hv, ldata& datalst) const
// 誤差の2階偏導関数
{
  double d2q = 0.0;
  for (auto& cdatai : datalst)
  {
    double x = cdatai.fdata(0), y = cdatai.fdata(1);
    if (x != BLANKf && y != BLANKf)
    {
      double qd2 = funcd2(x, hiidx, hjidx, hv);
      double d2ei = y * qd2 - funcd(x, hiidx, hv) * funcd(x, hjidx, hv) - func(x, hv) * qd2;
      d2q -= 2.0 * d2ei;
    }
  }
  return d2q;                   // 誤差2階導関数値を返す
}

double option::func(double x, const vec<double>& hv) const
// モデル関数
{
  double f = 0.0;
  switch (model())
  {
  case modelexp : f = hv(0) * exp(hv(1) * x); // 指数関数
    break;
  case modellog : f = logdis(x, hv(2), hv(1), hv(0)); // ロジスティック関数
    break;
  case modellog1 : f = logdis(x, 1.0, hv(1), hv(0)); // ロジスティック関数(M=1固定)
    break;
  case modellogc : f = ym() + logdis(x, hv(2), hv(1), hv(0) - ym()); // ロジスティック関数(最小値有)
    break;
  default : f = hv(0) + hv(1) * x; // 1次関数
  }
  return f;                     // 関数値を返す
}

double option::funcd(double x, int hidx, const vec<double>& hv) const
// モデル偏導関数
{
  constexpr double Delta = EPSdef; // 微小変動の割合(差分法用)
  vec<double> hdv(hv);
  double dh = 0.0, grad = 0.0;
  switch (model())
  {
  case modelexp :               // 指数関数
    switch (hidx)
    {
    case 0 : grad = exp(hv(1) * x); // hv(0)で偏微分
      break;
    case 1 : grad = x * exp(hv(1) * x); // hv(1)で偏微分
    }
    break;
  case modeletc :               // 1次関数
    switch (hidx)
    {
    case 0 : grad = 1.0;        // hv(0)で偏微分
      break;
    case 1 : grad = x;          // hv(1)で偏微分
    }
    break;
  default : dh = hv(hidx) * Delta; // 差分法
    hdv.assign(hidx, hdv(hidx) + dh);
    grad = (sgn(dh) == 0) ? 0.0 : (func(x, hdv) - func(x, hv)) / dh;
  }
  return grad;                  // 偏導関数値を返す
}

double option::funcd2(double x, int hiidx, int hjidx, const vec<double>& hv) const
// モデル2階偏導関数
{
  constexpr double Delta = EPSdef; // 微小変動の割合(差分法用)
  vec hijv(hv), hjv(hv);
  double fhi = 0.0, fhj = 0.0, dhi = 0.0, dhj = 0.0, grad2 = 0.0;
  switch (model())
  {
  case modelexp :               // 指数関数
    switch (hiidx)
    {
    case 0 :
      switch (hjidx)
      {
      case 0 : grad2 = 0.0;
        break;
      case 1 : grad2 = x * exp(hv(1) * x);
      }
      break;
    case 1 :
      switch (hjidx)
      {
      case 0 : grad2 = x * exp(hv(1) * x);
        break;
      case 1 : grad2 = x * x * hv(0) * exp(hv(1) * x);
      }
    }
    break;
  case modeletc : grad2 = 0.0;  // 1次関数
    break;
  default : dhi = hv(hiidx) * Delta; // 差分法
    dhj = hv(hjidx) * Delta;
    hijv.assign(hiidx, hv(hiidx) + dhi);
    fhi = func(x, hijv);
    hjv.assign(hjidx, hv(hjidx) + dhj);
    fhj = func(x, hjv);
    hijv.assign(hjidx, hijv(hjidx) + dhj);
    grad2 = (sgn(dhi * dhj) == 0) ? 0.0 : (func(x, hijv) - fhi - fhj + func(x, hv)) / dhi / dhj;
  }
  return grad2;                 // 2階偏導関数値を返す
}

double option::absmin(const vec<double>& hv) const
// 絶対値最小のパラメターの値を返す
{
  double hmin = FDmax;
  for (size_t i = 0;  i < hv.dim();  ++i)
  {
    double hi = fabs(hv(i));
    if (FDmin < hi && hi < hmin) hmin = hi;
  }
  return (hmin == FDmax) ? 1.0 : hmin;
}

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_nr;                     // 反復計算回数
  pairdbl xrange;               // xの最小値と最大値
  pairdbl yrange;               // yの最小値と最大値
  statcor<double> statcxy;      // xとyの相関基礎統計量
  vec<double> v_h0v;            // パラメータ初期値
  vec<double> v_hv;             // パラメータ推測値
  double v_qerr;                // 回帰残差
  restype v_resvar;             // 計算結果の種類
  restype init_var(option& opt, ldata& datalst); // データを読み込み変数を初期化する
  void init_h0(statcor<double>& statcxln, statcor<double>& statcxly, option& opt, ldata& datalst);
  // パラメータの初期値を設定する
  restype init_log(option& opt, ldata& datalst); // 微分回帰分析による初期値の設定
  int cal_smat(mat<double>& smat, double ym, ldata& datalst) const;
  // 微分回帰分析用に積和行列を計算する
  restype cal_simplex(int& nr, option& opt, ldata& datalst);
  // Simplex法によってパラメータを求める
  void init_simplex(mat<double>& simplexmat, const vec<double>& hv, option& opt) const;
  // シンプレックスの初期設定
  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, option& opt, ldata& datalst) const;
  // fiの計算と終了基準のチェック
  restype check_hv(const vec<double>& hv, option& opt) const;
  // パラメータの発散チェック
  restype cal_newton(int& nr, option& opt, ldata& datalst);
  // Newton-Raphson法によってパラメータを求める
  void cal_gh(vec<double>& gradv, mat<double>& hessmat, const vec<double>& hv, option& opt, ldata& datalst) const;
  // 傾斜ベクトルとHessの行列の計算
  restype check_q(double qnew, double qold, const vec<double>& hv, int nr, option& opt) const;
  // Q(h)の計算と終了基準のチェック
public:
  varclass(int nh0) : v_nr(0), xrange(FDmax, -FDmax), yrange(FDmax, -FDmax), statcxy(), v_h0v(nh0), v_hv(nh0), v_qerr(0.0), v_resvar(resnon)
  // デフォールトコンストラクタ
  {}
  ~varclass()                   // デストラクタ
  {}
  void cal(option& opt, ldata& datalst);
  // 非線形最小2乗法によってパラメータを求める
  size_t nh() const             // パラメータ数を返す
  { return v_h0v.dim(); }
  int nr() const                // 反復計算回数を返す
  { return v_nr; }
  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(); }
  double xysyy() const          // yの平方和を返す
  { return statcxy.syy(); }
  const vec<double>& h0v() const // パラメーター初期値を返す
  { return v_h0v; }
  double h0v(size_t idx) const  // 特定のパラメーター初期値を返す
  { return (idx < nh()) ? v_h0v(idx) : 0.0; }
  const vec<double>& hv() const // パラメーターを返す
  { return v_hv; }
  double hv(size_t idx) const   // 特定のパラメーターを返す
  { return (idx < nh()) ? v_hv(idx) : 0.0; }
  double qerr() const           // 回帰残差を返す
  { return v_qerr; }
  restype resvar() const        // 計算結果の種類を返す
  { return v_resvar; }
};

void varclass::cal(option& opt, ldata& datalst)
// 非線形最小2乗法によってパラメータを求める
{
  v_nr = 0;  v_resvar = resnon;
  if (init_var(opt, datalst) == reserr) // 初期値の設定失敗
  {
    v_resvar = reserr;
    return;
  }
  if (opt.meth() == methetc || opt.meth() == methsimplex)
  {
    v_resvar = cal_simplex(v_nr, opt, datalst);
    if (v_resvar == reserr)     // 収束失敗
    {
      v_qerr = opt.fq(v_hv, datalst);
      return;
    }
  }
  if (opt.meth() == methetc || opt.meth() == methnewton)
  {
    int nr = 0;
    v_resvar = cal_newton(nr, opt, datalst);
    v_nr += nr;
    if (v_resvar == reserr)     // 収束失敗
    {
      v_qerr = opt.fq(v_hv, datalst);
      return;
    }
  }
  v_qerr = opt.fq(v_hv, datalst);
}

restype varclass::init_var(option& opt, ldata& datalst)
// データを読み込み変数を初期化する
{
  statcor<double> statcxln;     // xとln(y)の相関基礎統計量
  statcor<double> statcxly;     // xとlogit(y)の相関基礎統計量
  for (auto& cdatai : datalst)
  {
    double x = cdatai.fdata(0), y = cdatai.fdata(1);
    if (x != BLANKf && y != 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;
      statcxy.sigma(x, y);
      if (0.0 < y) statcxln.sigma(x, log(y));
      if (0.0 < y && y < 1.0) statcxly.sigma(x, log(y / (1.0 - y)));
    }
  }
  if (statcxy.n() < opt.nh()) return reserr; // データ不足！
  init_h0(statcxln, statcxly, opt, datalst); // パラメータ初期値設定
  return resok;                              // 初期値設定完了
}

void varclass::init_h0(statcor<double>& statcxln, statcor<double>& statcxly, option& opt, ldata& datalst)
// パラメータの初期値を設定する
{
  double y0 = (yrange.first > 0.0) ? yrange.first * 0.5 : 0.0;
  double m = yrange.second;
  double b = (statcxy.sxx() > 0.0) ? statcxy.sxy() / statcxy.sxx() : 0.0;
  double c = (sgn(m) == 0) ? 0.0 : 8.0 * b / m; // ロジスティック関数の回帰係数
  double a = statcxy.my() - b * statcxy.mx();
  double bln = (statcxln.sxx() > 0.0) ? statcxln.sxy() / statcxln.sxx() : 0.0;
  double aln = statcxln.my() - bln * statcxln.mx();
  double bly = (statcxly.sxx() > 0.0) ? statcxly.sxy() / statcxly.sxx() : 0.0;
  double aly = bly * statcxly.mx() - statcxly.my();
  switch (opt.model())
  {
  case modellog :               // ロジスティック関数
    if (init_log(opt, datalst) == reserr) // 微分回帰分析による初期値設定失敗
    {
      v_h0v.assign(0, y0);
      v_h0v.assign(1, c);
      v_h0v.assign(2, m);
    }
    break;
  case modellog1 :              // ロジスティック関数(M=1固定)
    if (statcxly.n() > 1)       // ロジット回帰分析による初期値設定
    {
      y0 = 1.0 / (0.5 * exp(aly) + 1.0);
      c = bly;
      v_h0v.assign(0, y0);
      v_h0v.assign(1, bly);
    }
    else
    {
      y0 = 0.5 * a;
      c = 8.0 * b;
      v_h0v.assign(0, y0);
      v_h0v.assign(1, c);
    }
    break;
  case modellogc :              // ロジスティック関数(最小値有)
    if (init_log(opt, datalst) == reserr) // 微分回帰分析による初期値設定失敗
    {
      v_h0v.assign(0, y0);
      v_h0v.assign(1, c);
      v_h0v.assign(2, m);
      if (opt.ym() == BLANKf) opt.ym(0.5 * y0); // 最小値自動設定
    }
    break;
  case modelexp : v_h0v.assign(0, exp(aln)); // 指数関数
    v_h0v.assign(1, bln);
    break;
  default : v_h0v.assign(0, a);  v_h0v.assign(1, b); // その他の関数
  }
  if (opt.nh0() > 0)            // 初期値が指定されている場合
  {
    for (size_t i = 0; i < opt.nh0(); ++i) v_h0v.assign(i, opt.h0(i));
  }
  v_hv = v_h0v;
}

restype varclass::init_log(option& opt, ldata& datalst)
// 微分回帰分析による初期値の設定
{
  int p = Pdef;
  bool initym = false;
  double ym = 0.0;
  if (opt.model() == modellogc) // 最小値があるロジスティック関数
  {
    if (opt.ym() == BLANKf) { initym = true;  p = Pdef + 1; } // 最小値自動設定
    else ym = opt.ym();
  }
  size_t dim = p + 1;
  mat<double> smat(dim);        // 積和行列
  int nxy = cal_smat(smat, ym, datalst);
  if (nxy - 2 < p) return reserr; // 微分回帰分析できず
  size_t dyidx = p;             // dyのインデックス
  double b0 = 0.0, b1 = 0.0, b2 = 0.0, c = 0.0, m = 0.0;
  if (initym)
  {
    b0 = smat(0, dyidx);  b1 = smat(1, dyidx);  b2 = smat(2, dyidx);
    if (sgn(b1 * b2) >= 0) return reserr; // 2次回帰分析結果不適切
    double tmp = b1 * b1 - 4.0 * b2 * b0;
    tmp = (sgn(tmp) > 0) ? sqrt(tmp) : 0.0;
    ym = (-b1 + tmp) / (2.0 * b2);
    double ymax = (-b1 - tmp) / (2.0 * b2);
    if (sgndiff(ym, ymax) < 0)
    { m = ymax - ym; }
    else
    {
      m = ym - ymax;
      ym = ymax;
    }
    c = -b2 * m;  opt.ym(ym);
  }
  else
  {
    b1 = smat(0, dyidx);  b2 = smat(1, dyidx);
    if (sgn(b1 * b2) >= 0) return reserr; // 2次回帰分析結果不適切
    c = b1;  m = -b1 / b2;
  }
  statcor<double> statcxly;     // xとlogit(y)の相関基礎統計量
  for (auto& cdatai : datalst)
  {
    double x = cdatai.fdata(0), y = cdatai.fdata(1);
    if (x != BLANKf && y != BLANKf)
    {
      double yi = y - ym;
      if (0.0 < yi && yi < m) statcxly.sigma(x, log(yi / (m - yi)));
    }
  }
  if (statcxly.n() < 1) return reserr; // ロジットなし
  double ly0 = statcxly.my() - c * statcxly.mx();
  double y0 = ym + m / (1.0 + exp(-ly0));
  v_h0v.assign(0, y0);  v_h0v.assign(1, c);  v_h0v.assign(2, m);
  return resok;                 // 初期値設定完了
}

int varclass::cal_smat(mat<double>& smat, double ym, ldata& datalst) const
// 微分回帰分析用に積和行列を計算する
{
  size_t dim = smat.r();
  size_t p = smat.r() - 1;      // 2次回帰式のパラメータ数
  vecdbl fdatavec(dim);
  size_t offset = 0;
  if (dim == Pdef + 2) { fdatavec[0] = 1.0;  offset = 1; }
  double x1 = 0.0, x2 = 0.0, y1 = 0.0, y2 = 0.0;
  size_t nxy = 0;
  for (auto& cdatai : datalst)
  {
    double x = cdatai.fdata(0), y = cdatai.fdata(1);
    if (x != BLANKf && y != BLANKf)
    {
      if (nxy < 1)
      { x1 = x;  y1 = y;  nxy = 1; }
      else if (nxy == 1)
      { x2 = x;  y2 = y;  nxy = 2; }
      else
      {
        double dx1 = x2 - x1, dx2 = x - x2;
        double w1 = (fabs(dx1) < FDeps) ? 0.0 : 1.0 / dx1;
        double w2 = (fabs(dx2) < FDeps) ? 0.0 : 1.0 / dx2;
        double b1 = (y2 - y1) * w1, b2 = (y - y2) * w2;
        w1 = fabs(w1);  w2 = fabs(w2);
        double dy = (fabs(w1 + w2) < FDeps) ? 0.0 : (b1 * w1 + b2 * w2) / (w1 + w2);
        double yi = y2 - ym;    // 定数項を引く
        double fd = yi;
        for (size_t i = 0; i < Pdef; ++i)
        {
          fdatavec[offset + i] = fd;
          fd *= yi;
        }
        fdatavec[p] = dy;
        for (size_t i = 0; i < dim; ++i)
        {
          fd = fdatavec[i];
          for (size_t j = i; j < dim; ++j) smat.add(i, j, fd * fdatavec[j]);
        }
        x1 = x2;  y1 = y2;
        x2 = x;  y2 = y;  ++nxy;
      }
    }
  }
  for (size_t i = 0; i < p; ++i)
  {
    for (size_t j = i + 1; j <= p; ++j) smat.assign(j, i, smat(i, j));
  }
  for (size_t i = 0; i < p; ++i) smat.sweep(i);
  return nxy;               // 例数を返す
}

restype varclass::cal_simplex(int& nr, option& opt, ldata& datalst)
// Simplex法によってパラメータを求める
{
  int p = nh();                     // パラメーター数
  mat<double> simplexmat(p, p + 1); // シンプレックスを行列で表す
  vec<double> hlv(v_hv), hmv(v_hv), hhv(v_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, v_hv, opt);
  while ((ressimp = check_simplex(simplexmat, hlv, hmv, hhv, fl, fs, fh, lidx, hidx, nr, opt, datalst)) == resnon)
  {
    vec<double> hrv(hmv);
    hrv *= 2.0;
    hrv -= hhv;
    double fr = opt.fq(hrv, datalst);
    if (fr < fl)                // 折り返し成功
    {
      vec<double> hev(hrv);
      hev *= 2.0;
      hev -= hmv;
      double fe = opt.fq(hev, datalst);
      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 hcv(hmv + hhv);
      hcv *= 0.5;
      double fc = opt.fq(hcv, datalst);
      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);
          }
        }
      }
    }
    ++nr;
  }
  if (ressimp == reserr || fabs(fl) > FLT_MAX) return reserr; // 発散
  v_hv = hlv;
  return ressimp;               // 収束または反復回数最大値
}

void varclass::init_simplex(mat<double>& simplexmat, const vec<double>& hv, option& opt) const
// シンプレックスの初期設定
{
  constexpr double Sqrt2 = sqrt(2.0); // コンパイル時定数
  constexpr double Rs = 0.1;          // シンプレックスの1辺の相対割合
  int p = hv.dim();
  double s = opt.absmin(hv) * 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, option& opt, ldata& datalst) const
// fiの計算と終了基準のチェック
{
  int p = hlv.dim();
  int eidx = p - 1;
  vecdbl fvec(p + 1, 0.0);
  fl = FDmax;  fh = -FDmax;
  for (int j = 0; j <= p; ++j)
  {
    vec<double> hjv(simplexmat.sub(0, j, eidx, j));
    double fj = opt.fq(hjv, datalst);
    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 (check_hv(hlv, opt) == reserr) return reserr; // 発散
  double ferr = fabs(fh - fl);
  if (ferr < fabs(fl * opt.eps())) return resok; // 収束
  if (opt.rmax() > 0 && nr >= opt.rmax()) 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;                // 収束せず
}

restype varclass::check_hv(const vec<double>& hv, option& opt) const
// パラメータの発散チェック
{
  restype reshv = resnon;
  for (size_t i = 0; reshv == resnon && i < hv.dim(); ++i)
  {
    if (fabs(hv(i)) > FDmax || isnan(hv(i))) reshv = reserr;
  }
  if (reshv != resnon) return reshv; // 発散
  if (opt.model() == modellog || opt.model() == modellog1 || opt.model() == modellogc)
  {
    if (opt.model() == modellog1)
    {
      if (sgn(hv(0)) <= 0 || 1.0 <= hv(0)) reshv = reserr;
    }
    else if (opt.model() == modellog)
    {
      if (hv(0) >= hv(2) || hv(0) > yrange.second) reshv = reserr;
    }
    else if (opt.model() == modellogc)
    {
      if (hv(0) <= opt.ym() || hv(0) - opt.ym() >= hv(2) || hv(0) > yrange.second) reshv = reserr;
    }
  }
  return reshv;                 // チェック結果を返す
}

restype varclass::cal_newton(int& nr, option& opt, ldata& datalst)
// Newton-Raphson法によってパラメータを求める
{
  int p = nh();                 // パラメーター数
  vec<double> gradv(p), htv(v_hv);
  mat<double> hessmat(p, p);
  double qnew = opt.fq(htv, datalst), qold = 0.0;
  restype resnew = resnon;
  nr = 0;
  do
  {
    ++nr;
    cal_gh(gradv, hessmat, htv, opt, datalst);
    hessmat.inv();
    htv -= hessmat * gradv;
    qold = qnew;
    qnew = opt.fq(htv, datalst);
  } while ((resnew = check_q(qnew, qold, htv, nr, opt)) == resnon);
  if (resnew == reserr) return resnew; // 発散
  v_hv = htv;
  return resnew;                // 収束または反復回数最大値
}

void varclass::cal_gh(vec<double>& gradv, mat<double>& hessmat, const vec<double>& hv, option& opt, ldata& datalst) const
// 傾斜ベクトルとHessの行列の計算
{
  for (size_t i = 0; i < hv.dim(); ++i)
  {
    gradv.assign(i, opt.fqd(i, hv, datalst));
    for (size_t j = i; j < hv.dim(); ++j)
    { hessmat.assign(i, j, opt.fqd2(i, j, hv, datalst)); }
  }
  for (size_t i = 0; i < hv.dim() - 1; ++i)
  {
    for (size_t j = i + 1; j < hv.dim(); ++j) hessmat.assign(j, i, hessmat(i, j));
  }
}

restype varclass::check_q(double qnew, double qold, const vec<double>& hv, int nr, option& opt) const
// Q(h)の計算と終了基準のチェック
{
  if (check_hv(hv, opt) == reserr || qnew > FDmax || isnan(qnew)) return reserr; // 発散
  double qerr = fabs(qold - qnew);
  if (qerr < fabs(qnew * opt.eps())) return resok; // 収束
  if (opt.rmax() > 0 && nr >= opt.rmax()) return resrmax; // 反復計算回数最大値
  return resnon;                // 収束せず
}

class nlraclass                 // メインクラス
{
private:
  string fname;                 // ファイル名
  idd id;                       // データ識別情報
  ldata datalst;                // データリスト
  void put_resid(ostream& fout); // 分析結果のタイトルとデータ識別情報を出力する
  void put_model(ostream& fout, option& opt, const varclass& var) const;
  // モデル関数の非線形最小2乗法計算結果の出力
  void put_nlatab(ostream& fout, option& opt, const varclass& var);
  // 非線形最小2乗法計の分散分析表出力
  void put_exp(ostream& fout, option& opt, const varclass& var) const;
  // 指数関数の結果出力
  void put_log(ostream& fout, option& opt, const varclass& var) const;
  // ロジスティック関数の結果出力
  void put_log1(ostream& fout, option& opt, const varclass& var) const;
  // ロジスティック関数(M=1に固定)の結果出力
  void put_logc(ostream& fout, option& opt, const varclass& var) const;
  // ロジスティック関数(定数項有)の結果出力
  void put_etc(ostream& fout, option& opt, const varclass& var) const;
  // 1次関数の結果出力
  void plot_ngp(option& opt, const varclass& var);
  // Ngraph用データ出力
  void plot_ngpid(ofstream& fout);
  // Ngraph用データのタイトルとデータ識別情報を出力する
  void plot_data(ofstream& fout); // Ngraph用に個々のデータを出力する
  void plot_func(ofstream& fout, option& opt, const varclass& var);
  // Ngraph用関数曲線の出力
public:
  nlraclass() : fname(), id(), datalst() // コンストラクタ
  {}
  ~nlraclass()                  // デストラクタ
  {}
  int get(const char* c_fname, option& opt); // データを読み込み保存・累積する
  void put(option& opt);                     // 計算結果の出力
};

void nlraclass::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 nlraclass::put_model(ostream& fout, option& opt, const varclass& var) const
// モデル関数の非線形最小2乗法計算結果の出力
{
  switch (opt.model())
  {
  case modelexp : put_exp(fout, opt, var);
    break;
  case modellog : put_log(fout, opt, var);
    break;
  case modellog1 : put_log1(fout, opt, var);
    break;
  case modellogc : put_logc(fout, opt, var);
    break;
  default : put_etc(fout, opt, var);
  }
  int n = var.xyn(), p = var.nh();
  if (n < p)                    // データ不足！
  {
    fout << "※関数計算不可能！" << endl;
    return;
  }
  double aic = (sgn(var.qerr()) > 0) ? (double)n * log(var.qerr()) + 2.0 * (double)p : 0.0;
  double sa = var.xysyy(), sl = var.qerr();
  if (sl < FDeps) sl = 0.0;
  double sc = sa - sl;
  if (sc < FDeps) sc = 0.0;
  double r2 = (sa > 0.0) ? sc / sa : 0.0;
  fout << "\n・計算方法:" << opt.methname() << "  反復回数=" << var.nr();
  if (var.resvar() == reserr)
  { fout << ":非線形最小2乗法発散！"; }
  else
  {
    if (var.resvar() == resrmax) fout << "(最大反復回数)";
    fout << "\n  残差平方和(Q)=" << var.qerr() << "  赤池の情報量基準(AIC)=" << aic << "\n";
    fout << "  重寄与率:R2乗=" << r2 << "  重相関係数:R=" << sqrt(r2);
  }
  fout << endl;
}

void nlraclass::put_nlatab(ostream& fout, option& opt, const varclass& var)
// 非線形最小2乗法計の分散分析表出力
{
  int n = var.xyn(), p = var.nh();
  double sa = var.xysyy(), sl = var.qerr();
  if (sl < FDeps) sl = 0.0;
  int dfa = n - 1, dfc = p - 1;
  double sc = sa - sl;
  if (sc < FDeps) sc = 0.0;
  int dfl = dfa - dfc;
  if (opt.res("atab"))
  {
    fout << "\n";
    putatab(fout, "回帰", sc, dfc, sl, dfl, sa, dfa);
  }
  if (opt.res("def") == false && opt.res("matab") == false) return;
  vec htv(var.hv());
  switch (opt.model())
  {
  case modellog : htv.assign(0, var.hv(2) * 0.5);  htv.assign(1, 0.0);
    sa = opt.fq(htv, datalst);
    break;
  case modellog1 : htv.assign(0, 0.5);  htv.assign(1, 0.0);
    sa = opt.fq(htv, datalst);
    break;
  case modellogc : htv.assign(0, var.hv(2) * 0.5 + opt.ym());  htv.assign(1, 0.0);
    sa = opt.fq(htv, datalst);
    break;
  default : return;
  }
  sc = sa - sl;
  if (sc < FDeps) sc = 0.0;
  double r2 = (sa > 0.0) ? sc / sa : 0.0;
  fout << "\n・全体の平方和をモデル固有のものにした時\n";
  fout << "  重寄与率:R2乗=" << r2 << "  重相関係数:R=" << sqrt(r2) << "\n\n";
  putatab(fout, "回帰", sc, dfc, sl, dfl, sa, dfa);
}

void nlraclass::put_exp(ostream& fout, option& opt, const varclass& var) const
// 指数関数の結果出力
{
  fout << "\n・モデル:" << opt.modelname() << "\n";
  fout << "  モデル関数:y=A・exp(c・x)\n";
  fout << "  パラメーター初期値:A=" << var.h0v(0) << "  c=" << var.h0v(1) << "\n";
  fout << "  パラメーター推定値:A=" << var.hv(0) << "  c=" << var.hv(1) << endl;
}

void nlraclass::put_log(ostream& fout, option& opt, const varclass& var) const
// ロジスティック関数の結果出力
{
  double y00 = var.h0v(0), c0 = var.h0v(1), m0 = var.h0v(2);
  double y0 = var.hv(0), c = var.hv(1), m = var.hv(2);
  double ly00 = (sgn(m0 - y00) == 0) ? 0.0 : log(y00 / (m0 - y00));
  double ly0 = (sgn(m - y0) == 0) ? 0.0 : log(y0 / (m - y0));
  fout << "\n・モデル:" << opt.modelname() << "\n";
  fout << "  モデル関数:y=M/[1+{(M-y0)/y0}exp(-c・x)]=M/{1+exp(-ly0-c・x)}\n";
  fout << "  パラメーター初期値:y0=" << y00 << "  c=" << c0 << "  M=" << m0 << "  ly0=" << ly00 << "\n";
  fout << "  パラメーター推定値:y0=" << y0 << "  c=" << c << "  M=" << m << "  ly0=" << ly0 << endl;
}

void nlraclass::put_log1(ostream& fout, option& opt, const varclass& var) const
// ロジスティック関数(M=1固定)の結果出力
{
  double y00 = var.h0v(0), c0 = var.h0v(1);
  double y0 = var.hv(0), c = var.hv(1);
  double ly00 = (sgn(1.0 - y00) == 0) ? 0.0 : log(y00 / (1.0 - y00));
  double ly0 = (sgn(1.0 - y0) == 0) ? 0.0 : log(y0 / (1.0 - y0));
  fout << "\n・モデル:" << opt.modelname() << "\n";
  fout << "  モデル関数:y=1/[1+{(1-y0)/y0}exp(-c・x)]=1/{1+exp(-ly0-c・x)}\n";
  fout << "  パラメーター初期値:y0=" << y00 << "  c=" << c0 << "  ly0=" << ly00 << "\n";
  fout << "  パラメーター推定値:y0=" << y0 << "  c=" << c << "  ly0=" << ly0 << endl;
}

void nlraclass::put_logc(ostream& fout, option& opt, const varclass& var) const
// ロジスティック関数(最小値有)の結果出力
{
  double y00 = var.h0v(0), c0 = var.h0v(1), m0 = var.h0v(2);
  double y0 = var.hv(0), c = var.hv(1), m = var.hv(2);
  double ly00 = (sgn(m0 - (y00 - opt.ym())) == 0) ? 0.0 : log((y00 - opt.ym()) / (m0 - (y00 - opt.ym())));
  double ly0 = (sgn(m - (y0 - opt.ym())) == 0) ? 0.0 : log((y0 - opt.ym()) / (m - (y0 - opt.ym())));
  fout << "\n・モデル:" << opt.modelname() << "\n";
  fout << "  モデル関数:y=ym+M/[1+[{M-(y0-ym)}/(y0-ym)]exp(-c・x)]\n";
  fout << "              =ym+M/{1+exp(-ly0-c・x)}\n";
  fout << "  パラメーター初期値:y0=" << y00 << "  c=" << c0 << "  m=" << m0 << "\n";
  fout << "                     ym=" << opt.ym() << "  ly0=" << ly00 << "\n";
  fout << "  パラメーター推定値:y0=" << y0 << "  c=" << c << "  m=" << m << "\n";
  fout << "                     ym=" << opt.ym() << "  ly0=" << ly0 << endl;
}

void nlraclass::put_etc(ostream& fout, option& opt, const varclass& var) const
// 1次関数の結果出力
{
  fout << "\n・モデル:" << opt.modelname() << "\n";
  fout << "  モデル関数:y=a + b・x\n";
  fout << "  パラメーター初期値:a=" << var.h0v(0) << "  b=" << var.h0v(1) << "\n";
  fout << "  パラメーター推定値:a=" << var.hv(0) << "  b=" << var.hv(1) << endl;
}

void nlraclass::plot_ngp(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);
    plot_data(fout);
    delml = true;
  }
  if (opt.ngp("func"))
  {
    if (delml) fout << NGPdelm << "\n";
    else fout << "-s1 -dline\n";
    plot_ngpid(fout);
    plot_func(fout, opt, var);
  }
}

void nlraclass::plot_ngpid(ofstream& fout)
// 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";
}

void nlraclass::plot_data(ofstream& fout)
// Ngraph用に個々のデータを出力する
{
  fout << "%生データ\n";
  for (auto& cdatai : datalst)
  {
    double x = cdatai.fdata(0), y = cdatai.fdata(1);
    if (x != BLANKf && y != BLANKf) fout << x << " " << y << endl;
  }
}

void nlraclass::plot_func(ofstream& fout, option& opt, const varclass& var)
// Ngraph用関数曲線の出力
{
  fout << "%モデル:" << opt.modelname() << "\n";
  fout << "%計算方法:" << opt.methname() << "  反復回数=" << var.nr();
  if (var.resvar() == resrmax) fout << "(最大反復回数)";
  fout << "\n";
  double xmin = var.xmin(), xmax = var.xmax(), xmg = 0.1;
  if (opt.xmin() != BLANKf && opt.xmax() != BLANKf)
  { xmin = opt.xmin();  xmax = opt.xmax();  xmg = 0.0; }
  double xrange = xmax - xmin;
  xmin -= xrange * xmg;  xmax += xrange * xmg;
  double dx = xrange / (double)RDdef;
  double x = xmin;
  for (int i = 0; i <= RDdef; ++i, x += dx)
  { fout << x << " " << opt.func(x, var.hv()) << endl; }
}

int nlraclass::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 nlraclass::put(option& opt) // 計算結果の出力
{
  datalst.sort(data_less(0));   // 説明変数順にソートする
  varclass var(opt.nh());       // 変数オブジェクト
  var.cal(opt, datalst);
  constexpr int LENline = 69;   // 区切り線の長さ
  if (opt.res("res"))           // 計算結果出力
  {
    put_resid(cout);
    if (opt.res("stat"))
    {
      cout << "\n・変数の基礎統計量\n";
      putline(cout, LENline);
      int n = var.xyn();
      cout << "x:";  putmsd(cout, n, var.xymx(), var.xyvxx());
      cout << "y:";  putmsd(cout, n, var.xymy(), var.xyvyy());
      putline(cout, LENline);
    }
    if (var.xyn() < 1)
    {
      cout << "※有効なデータなし！\n\n" << endl;
      return;
    }
    if (opt.res("model")) put_model(cout, opt, var);
    if (var.resvar() != reserr) put_nlatab(cout, opt, var);
    cout << "\n" << endl;
  }
  if (var.resvar() != reserr && opt.ngpname() != "") plot_ngp(opt, var);
}

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