// mlog.C Version 7.44  Copyleft (c) 杉本典夫
// Time-stamp: <2025-10-31 17:02:33 tac>
//
// 累積ロジットモデル(比例オッズモデル)による順序ロジスティック回帰分析。
// ロジスティック回帰分析における群項目を順序のある多分類に拡張した手法。
//
// 使用ライブラリ:libdans.so, libfun.so
// コンパイル法 :g++ -O -s -W mlog mlog.C
//               -I/usr/local/include -L/usr/local/lib -ldans -lfun

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

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

const vecstr USAGE =            // プログラム使用法表示用
{
  "mlog  Version 7.44           Copyleft(Vo\\oV) 2025 杉本典夫",
  "書式:mlog [<オプション> …]  [<データファイル> …]",
  "機能:最初の項目を群項目(順序分類データ)、以後の項目を説明変数Xとして",
  "     累積ロジットモデルによる順序ロジスティック回帰分析を行います。",
  "<オプション>",
  "  -b=<β0の初期値>[,<β1の初期値>,…]:偏回帰係数βの初期値を指定します。",
  "    オプション省略時:判別分析に基づいた初期値を自動計算",
  "  -cc[=<信頼係数>]:信頼係数を指定します。(0≦信頼係数<1)",
  "    -cc:-cc=0.95(オプション省略時)",
  "  -dnsname[=<出力ファイル名>]:DANSデータを出力するファイル名を指定する。",
  "    オプション省略時:DANSデータ出力せず  -dnsname:-dnsname=mlog_dns.d",
  "  -eps[=<収束誤差>]:最尤法の収束誤差を指定します。(0<収束誤差)",
  "  -f=<ファイル>:指定したファイルからオプションを読み込みます。",
  "  -h[=<コマンド名>]:オンラインマニュアルを表示します。(-h:-h=man)",
  "  -ngp={log|score}[,…]:Ngraph用出力データを指定します。",
  "    log:ロジスティック曲線  score:ロジットスコア",
  "  -ngpname[=<出力ファイル名>]:Ngraph用データを出力するファイル名を指定します。",
  "    オプション省略時:Ngraph用データ出力せず  -ngpname:-ngpname=mlog_ngp.dat",
  "  -res=[-][{def|all|ltab|olra|stat|v}][,…]:解析結果の出力項目を指定します。",
  "    -:-を付けた項目を出力しない(-res=-:全ての項目を出力しない)",
  "    def:オプション省略時と同じ項目を出力する  all:全ての項目を出力する",
  "    ltab:尤度分析表  olra:順序ロジスティック回帰分析  stat:基礎統計量",
  "    v:偏回帰係数初期値",
  "  -rmax=<最大回数>:反復計算の最大回数を指定します。(0≦最大回数)",
  "  -xrange=<最小値>-<最小値>:X軸の最小値と最大値を指定します。"
};

const string Title = "=== 順序ロジスティック回帰分析 ==="; // 表題
const string DNSNAMEdef = "mlog_dns.d"; // 省略時DANSファイル名
const string NGPNAMEdef = "mlog_ngp.dat"; // 省略時Ngraphファイル名
constexpr double EPSdef = 1.0e-5; // 収束誤差省略時値
constexpr int RMAXdef = 200;      // 繰り返し系算数省略時値

class option                    // オプションクラス
{
private:
  string v_path;                // プログラムのパス
  string v_pname;               // プログラム名
  lstr blst;                    // 説明変数の初期値リスト
  vecdbl v_bvec;                // 説明変数の初期値ベクトル
  double v_cc;                  // 信頼係数
  string v_dnsname;             // DANSファイル名
  double v_eps;                 // 最尤法の収束誤差
  string v_ngpname;             // Ngraphファイル名
  lstr ngplst;                  // Ngraph出力リスト
  mapsb ngpmap;                 // Ngraph出力の有無マップ
  lstr reslst;                  // 分析結果出力リスト
  mapsb resmap;                 // 分析結果出力の有無マップ
  int v_rmax;                   // 反復計算の最大回数
  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; }
  const vecdbl& bvec() const    // 説明変数の初期値ベクトルを返す
  { return v_bvec; }
  double b(size_t idx) const    // 説明変数の初期値を返す
  { return (idx < v_bvec.size()) ? v_bvec[idx] : 0.0; }
  double cc() const             // 信頼係数を返す
  { return v_cc; }
  const string& dnsname() const // DANSデータファイル名を返す
  { return v_dnsname; }
  double eps() const            // 最尤法の収束誤差を返す
  { return v_eps; }
  const string& ngpname() const // Ngraphデータファイル名を返す
  { return v_ngpname; }
  bool ngp(const string& code); // Ngraph出力フラッグを返す
  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; }
  void setdef()                 // 省略時オプション指定を設定する
  { getoptdef(*this); }
  void set(const string& optstr, char devc = OPTdev);
  // オプション指定によってオプションを設定する
  int check(int p);             // オプションのチェックとセット
  int set_ngp();                // -ngpオプションのセット
  int set_res();                // -resオプションのセット
};

option::option(const char* c_pname) : v_path(), v_pname(), blst(), v_bvec(), v_cc(CCdef), v_dnsname(), v_eps(EPSdef), v_ngpname(), ngplst(), ngpmap(), reslst(), resmap(), v_rmax(RMAXdef), 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] == "b")            // 偏回帰係数初期値の指定
  {
    if (strsplit(blst, opt[1], CHRdelm) < 1) err = badopt(pname(), optstr, USAGE);
  }
  else 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] == "dnsname") // DANSファイル名指定
  {
    if (opt[1].empty()) v_dnsname = DNSNAMEdef;
    else v_dnsname = opt[1];
    if (v_dnsname.size() < 3) err = badopt(pname(), optstr, USAGE);
  }
  else 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] == "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] == "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(int p)        // オプションのチェックとセット
{
  if (set_ngp() < 0) return -2; // -ngpオプションのセット
  if (set_res() < 0) return -3; // -resオプションのセット
  if (blst.empty()) return 0;   // 偏回帰係数初期値の指定なし
  if ((int)(blst.size()) > p) return disperr(pname(), "偏回帰係数初期値の値が多すぎます！");
  v_bvec.resize(p, BLANKf);
  int i = 0;
  for (const string& bstr : blst)
  {
    double fdata = strtodf(bstr);
    if (fdata == BLANKf) return disperr(pname(), "初期値指定が間違っています→", bstr);
    v_bvec[i] = fdata;
    ++i;
  }
  return 0;                     // エラーなし
}

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

class gstat                     // 群別統計量クラス
{
private:
  mat<double> v_smat;           // 積和行列
  int v_n;                      // 例数
  double v_my;                  // yの平均値
public:
  gstat() : v_smat(), v_n(0), v_my(BLANKf)
  // デフォールトコンストラクタ
  {}
  explicit gstat(size_t ni) : v_smat(ni + 1), v_n(0), v_my(BLANKf)
  // コンストラクタ
  {}
  gstat(const gstat& gstat0) : v_smat(gstat0.v_smat), v_n(gstat0.v_n), v_my(gstat0.v_my)
  // コピーコンストラクタ
  {}
  virtual ~gstat()              // デストラクタ
  {}
  const gstat& operator=(const gstat& gstat0); // コピー代入演算子
  const gstat& operator+=(const gstat& gstat0); // 加算代入演算子
  const mat<double>& smat() const               // 積和行列を返す
  { return v_smat; }
  double smat(int ridx, int cidx) const // 積和行列の要素を返す
  { return v_smat(ridx, cidx); }
  int n() const                 // 例数を返す
  { return v_n; }
  double my() const             // yの平均値を返す
  { return v_my; }
  double syy() const            // yの平方和を返す
  { return v_smat(v_smat.r() - 1, v_smat.c() - 1); }
  double sr() const             // 回帰残差を返す
  { return syy(); }
  void sigma(const vecdbl& fdata); // データを累積する
  int setsmat(vecdbl& mevec, vecdbl& vrvec);
  // 積和行列・平均値ベクトル・分散ベクトルを設定して例数を返す
};

const gstat& gstat::operator=(const gstat& gstat0) // コピー代入演算子
{
  if (&gstat0 == this) return *this; // 自分自身への代入を防ぐ
  v_smat = gstat0.v_smat;  v_n = gstat0.v_n;  v_my = gstat0.v_my;
  return *this;
}

const gstat& gstat::operator+=(const gstat& gstat0) // 加算代入演算子
{
  v_smat += gstat0.v_smat;  v_n += gstat0.v_n;
  return *this;
}

void gstat::sigma(const vecdbl& fdatavec) // データを累積する
{
  for (size_t i = 0; i < fdatavec.size(); ++i)
  {
    double x = fdatavec[i];
    for (size_t j = i; j < fdatavec.size(); ++j) v_smat.add(i, j, x * fdatavec[j]);
  }
}

int gstat::setsmat(vecdbl& mevec, vecdbl& vrvec)
// 積和行列・平均値ベクトル・分散ベクトルを設定して例数を返す
{
  v_n = (int)(v_smat(0, 0));
  for (size_t i = 0; i < v_smat.r() - 1; ++i)
  {
    for (size_t j = i + 1; j < v_smat.c(); ++j) v_smat.assign(j, i, v_smat(i, j));
  }
  v_smat.sweep(0);
  for (size_t i = 0; i < v_smat.r(); ++i)
  {
    mevec[i] = v_smat(0, i);
    vrvec[i] = (v_n > 1) ? v_smat(i, i) / (double)(v_n - 1) : 0.0;
  }
  v_my = v_smat(0, v_smat.r() -1);
  return v_n;                   // 例数を返す
}

using mapdg = map<double, gstat>;  // 群別統計量マップの別名定義
using mapdg_itr = mapdg::iterator; // 反復子の別名定義
using mapdg_citr = mapdg::const_iterator; // const反復子の別名定義

class varclass                 // 変数クラス
{
private:
  mapdg v_gsmap;                // 群統計量マップ
  gstat gstatall;               // 全体の統計量
  mat<double> smat0;            // 群0の判別分析用積和行列
  mat<double> smat1;            // 群1の判別分析用積和行列
  vecdbl v_mevec;               // 全体の平均値ベクトル
  vecdbl v_vrvec;               // 全体の分散ベクトル
  vec<double> b0v;              // 偏回帰係数初期値ベクトル
  vec<double> bv;               // 偏回帰係数ベクトル
  vec<double> vbv;              // 偏回帰係数の分散ベクトル
  double v_lb;                  // 対数尤度
  int v_ni;                     // 項目数
  int v_np;                     // パラメータ数
  int v_nt;                     // 総例数
  int v_nrep;                   // 反復計算数
  restype v_resvar;             // 処理結果の種類
  void set_initb(vec<double>& b0v, option& opt, const mapdg& gsmap, const gstat& gstatall, const mat<double>& smat0, const mat<double>& smat1) const;
  // 判別分析による偏回帰係数初期値のセット
  restype newton_mlog(int& nrep, vec<double>& bv, vec<double>& vbv, double& lb,  const mapdg& gsmap, option opt, ldata& datalst) const;
  // Newton法によるパラメーターの計算(最尤法)
  restype cal_lgh(vec<double>& gv, mat<double>& hmat, double& lb, const vec<double>& bv, const mapdg& gsmap, option opt, ldata& datalst) const;
  // 傾斜ベクトルgとHessian行列Hの計算
public:
  varclass(size_t ni) : v_gsmap(), gstatall(ni), smat0(ni + 1), smat1(ni + 1), v_mevec(ni + 1, 0.0), v_vrvec(ni + 1, 0.0), b0v(), bv(), vbv(), v_lb(0.0), v_ni(ni), v_np(0), v_nt(0), v_nrep(0), v_resvar(resnon)
  // コンストラクタ
  {}
  ~varclass()                   // デストラクタ
  {}
  const mapdg& gsmap() const    // 群統計量マップを返す
  { return v_gsmap; }
  const mat<double>& smat() const // 全体の積和行列を返す
  { return gstatall.smat(); }
  const vecdbl& mevec() const   // 全体の平均値ベクトルを返す
  { return v_mevec; }
  const vecdbl& vrvec() const   // 全体の分散ベクトルを返す
  { return v_vrvec; }
  double b0(int idx) const      // 偏回帰係数初期値を返す
  { return (0 <= idx && idx < np()) ? b0v(idx) : 0.0; }
  double b(int idx) const       // 偏回帰係数を返す
  { return (0 <= idx && idx < np()) ? bv(idx) : 0.0; }
  double vb(int idx) const      // 偏回帰係数の分散を返す
  { return (0 <= idx && idx < np()) ? vbv(idx) : 0.0; }
  double bst(int idx) const     // 標準偏回帰係数を返す
  { return (0 <= idx && idx < np()) ? b(idx) * sqrt(v_vrvec[idx - g() + 2]) : 0.0; }
  double lb() const             // 対数尤度を返す
  { return v_lb; }
  int ni() const                // 項目数を返す
  { return v_ni; }
  int np() const                // パラメータ数を返す
  { return v_np; }
  int g() const                 // 群数を返す
  { return (int)(v_gsmap.size()); }
  int nt() const                // 総例数を返す
  { return v_nt; }
  int nrep() const              // 反復計算回数を返す
  { return v_nrep; }
  restype resvar() const        // 処理結果の種類を返す
  { return v_resvar; }
  void cal(option& opt, ldata& datalst); // データを集計して各種統計量を求める
  void mlog(option& opt, ldata& datalst); // 順序ロジスティック回帰分析実行
  int search_index(double yg) const; // 目的変数データを検索し、インデックスを返す
  double zscore(int bidx, datac& cdata) const; // ロジットスコアを計算して返す
};

void varclass::cal(option& opt, ldata& datalst)
// データを集計して各種統計量を求める
{
  vecdbl fdatavec(ni() + 1, 0.0);
  fdatavec[0] = 1.0;
  for (auto& cdatai : datalst)
  {
    bool blank = false;
    for (int j = 1; blank == false && j < ni(); ++j)
    {
      fdatavec[j] = cdatai.fdata(j);
      if (fdatavec[j] == BLANKf) blank = true;
    }
    if (cdatai.fdata(0) == BLANKf || blank)
    { cdatai.assign(0, BLANKf); } // 後の集計のために欠損にする
    else
    {
      fdatavec[ni()] = cdatai.fdata(0);
      mapdg_itr gitr = v_gsmap.find(cdatai.fdata(0));
      if (gitr == v_gsmap.end())  // 新規群コード
      {
        gstat gstati(ni());
        gstati.sigma(fdatavec);
        v_gsmap.insert(make_pair(cdatai.fdata(0), gstati));
      }
      else                      // すでにある群
      { gitr->second.sigma(fdatavec); }
    }
  }
  if (v_gsmap.size() < 1)       // データなし
  {
    v_resvar = reserr;
    return;
  }
  int i = 1;
  for (auto& gstati : v_gsmap)
  {
    gstatall += gstati.second;
    if (i > 1) smat1 += gstati.second.smat(); // 2番目以後の単純積和行列を累積する
    gstati.second.setsmat(v_mevec, v_vrvec);     // 単純積和行列を積和行列にする
    if (i == 1) smat0 = gstati.second.smat(); // 最初の群の積和行列を保存する
    ++i;
  }
  v_nt = gstatall.setsmat(v_mevec, v_vrvec);
  for (size_t i = 0; i < smat1.r() - 1; ++i)
  {
    for (size_t j = i + 1; j < smat1.c(); ++j) smat1.assign(j, i, smat1(i, j));
  }
  smat1.sweep(0);               // 群1の単純積和行列→積和行列
  v_np = ni() + (int)(v_gsmap.size()) - 2;
  if (v_gsmap.size() < 2 || v_nt < v_np) v_resvar = reserr; // ロジスティック関数計算不可能
  else v_resvar = resok;
}

void varclass::mlog(option& opt, ldata& datalst)
// 順序ロジスティック回帰分析実行
{
  b0v.resize(np());  bv.resize(np());  vbv.resize(np());
  set_initb(b0v, opt, v_gsmap, gstatall, smat0, smat1);
  bv = b0v;
  v_resvar = newton_mlog(v_nrep, bv, vbv, v_lb, v_gsmap, opt, datalst);
}

void varclass::set_initb(vec<double>& b0v, option& opt, const mapdg& gsmap, const gstat& gstatall, const mat<double>& smat0, const mat<double>& smat1) const
// 判別分析による偏回帰係数初期値のセット
{
  mat<double> dmat(smat0 + smat1);
  dmat *= 1.0 / (double)(gstatall.n() - 2);
  dmat.assign(0, 0, 0.0);  dmat.assign(0, ni(), 0.0);
  dmat.assign(ni(), 0, 0.0);  dmat.assign(ni(), ni(), 0.0);
  for (int i = 1; i < ni(); ++i)
  {
    dmat.assign(0, i, (smat0(0, i) + smat1(0, i)) * 0.5);
    dmat.assign(i, 0, -dmat(0, i));
    dmat.assign(i, ni(), smat1(0, i) - smat0(0, i));
    dmat.assign(ni(), i, dmat(i, ni()));
  }
  for (int i = 1; i < ni(); ++i) dmat.sweep(i);
  mapdg_citr gitr = gsmap.begin();
  int n = gstatall.n(), nk = gitr->second.n();
  double cf = (n > nk && nk > 0) ? log((double)(n - nk) / (double)nk) : 0.0;
  b0v.assign(0, cf + dmat(0, ni()));
  ++gitr;
  for (int i = 1; i < g() - 1; ++i, ++gitr)
  {
    nk += gitr->second.n();
    double cfk = (n > nk && nk > 0) ? log((double)(n - nk) / (double)nk) : 0.0;
    b0v.assign(i, b0v(i - 1) + (cfk - cf));
    if (fabs(b0v(i)) < FDeps) b0v.assign(i, 0.0);
    cf = cfk;
  }
  for (int i = g() - 1; i < np(); ++i)
  {
    b0v.assign(i, dmat(i - g() + 2, ni()));
    if (fabs(b0v(i)) < FDeps) b0v.assign(i, 0.0);
  }
  if (opt.bvec().empty()) return; // 初期値指定なし
  for (int i = 0; i < np(); ++i)
  {
    if (opt.b(i) != BLANKf) b0v.assign(i, opt.b(i));
  }
}

restype varclass::newton_mlog(int& nrep, vec<double>& bv, vec<double>& vbv, double& lb, const mapdg& gsmap, option opt, ldata& datalst) const
// Newton法によるパラメーターの計算(最尤法)
{
  vec<double> gv(bv.dim());
  mat<double> hmat(bv.dim());
  double lold = 0.0;
  restype resnew = resnon;
  nrep = 0;
  while (resnew == resnon && cal_lgh(gv, hmat, lb, bv, gsmap, opt, datalst) != reserr)
  {
    ++nrep;
    hmat.inv();
    if (fabs(lb - lold) < opt.eps()) resnew = resok;
    else if (opt.rmax() > 0 && nrep >= opt.rmax()) resnew = resrmax;
    else { lold = lb;  bv -= hmat * gv; }
  }
  if (resnew == resnon) return reserr; // 発散
  for (size_t i = 0; i < bv.dim(); ++i) vbv.assign(i, -hmat(i, i));
  return resnew;                // Newton法成功
}

restype varclass::cal_lgh(vec<double>& gv, mat<double>& hmat, double& lb, const vec<double>& bv, const mapdg& gsmap, option opt, ldata& datalst) const
// 傾斜ベクトルgとHessian行列Hの計算
{
  constexpr double Emax = 308.0; // 指数最大値
  int p = (int)(bv.dim()), g = (int)(gsmap.size());
  vec<double> xiv(ni() - 1);
  gv = 0.0;  hmat = 0.0;  lb = 0.0;
  restype reslgh = resnon;
  ldata_itr ditr = datalst.begin();
  while (reslgh == resnon && ditr != datalst.end())
  {
    if (ditr->fdata(0) != BLANKf)
    {
      double bxi = 0.0;
      for (int j = 0; j < ni() - 1; ++j)
      {
        xiv.assign(j, ditr->fdata(j + 1));
        bxi += bv(g - 1 + j) * xiv(j);
      }
      int k = search_index(ditr->fdata(0));
      double bxik1 = (k == 0) ? 0.0 : bv(k - 1) + bxi;
      double bxik = (k == g - 1) ? 0.0 : bv(k) + bxi;
      if (fabs(bv(k)) > Emax || fabs(bxik1) > Emax || fabs(bxik) > Emax)
      { reslgh = reserr; }
      else
      {
        double eik1 = exp(bxik1), eik = exp(bxik), ei = 0.0, ei2 = 0.0;
        if (k == 0)
        {
          lb -= log(1.0 + eik);
          ei = eik / (1.0 + eik);  ei2 = eik / (1.0 + eik) / (1.0 + eik);
          gv.add(k, -ei);  hmat.add(k, k, -ei2);
          for (int j = g - 1; j < p; ++j) hmat.add(k, j, -xiv(j - g + 1) * ei2);
        }
        else if (k == g - 1)
        {
          lb += bxik1 - log(1.0 + eik1);
          ei = -1.0 / (1.0 + eik1);  ei2 = eik1 / (1.0 + eik1) / (1.0 + eik1);
          gv.add(k - 1, -ei);
          hmat.add(k - 1, k - 1, -ei2);
          for (int j = g - 1; j < p; ++j) hmat.add(k - 1, j, -xiv(j - g + 1) * ei2);
        }
        else
        {
          double d = eik1 / (1.0 + eik1) - eik / (1.0 + eik);
          lb += (d > 0.0) ? log(d) : 0.0;
          d = bv(k - 1) - bv(k);
          gv.add(k - 1, ((d != 0.0) ? 1.0 / (exp(d) - 1.0) : 0.0) + 1.0 / (1.0 + eik1));
          gv.add(k, ((d != 0.0) ? 1.0 / (exp(-d) - 1.0) : 0.0) + 1.0 / (1.0 + eik));
          d = exp(-bv(k - 1)) - exp(-bv(k));
          d = (d != 0.0) ? exp(-bv(k - 1) - bv(k)) / d / d : 0.0;
          hmat.add(k - 1, k - 1, -(d + eik1 / (1.0 + eik1) / (1.0 + eik1)));
          hmat.add(k, k, -(d + eik / (1.0 + eik) / (1.0 + eik)));
          hmat.add(k - 1, k, d);
          for (int j = g - 1; j < p; ++j)
          {
            hmat.add(k - 1, j, -xiv(j - g + 1) * eik1 / (1.0 + eik1) / (1.0 + eik1));
            hmat.add(k, j, -xiv(j - g + 1) * eik / (1.0 + eik) / (1.0 + eik));
          }
          ei = 1.0 - 1.0 / (1.0 + eik1) - 1.0 / (1.0 + eik);
          ei2 = eik1 / (1.0 + eik1) / (1.0 + eik1) + eik / (1.0 + eik) / (1.0 + eik);
        }
        for (int j = g - 1; j < p; ++j)
        {
          double xij = xiv(j - g + 1);
          gv.add(j, -xij * ei);
          for (int l = j; l < p; ++l) hmat.add(j, l, -xij * xiv(l - g + 1) * ei2);
        }
      }
    }
    ++ditr;
  }
  if (reslgh == reserr) return reslgh; // 計算エラー
  for (int i = 0; i < p - 1; ++i)
  { for (int j = i + 1; j < p; ++j) hmat.assign(j, i, hmat(i, j)); }
  return resok;
}

int varclass::search_index(double yg) const
// 目的変数データを検索し、インデックスを返す
{
  mapdg_citr gitr = v_gsmap.find(yg);
  if (gitr == v_gsmap.end()) return -1; // データなし
  return (int)(distance(v_gsmap.begin(), gitr)); // 先頭と要素の位置の差を返す
}

double varclass::zscore(int bidx, datac& cdata) const
// ロジットスコアを計算して返す
{
  double z = b(bidx);
  bool blank = false;
  for (int i = 1; blank == false && i < ni(); ++i)
  {
    double x = cdata.fdata(i);
    if (x == BLANKf) blank = true;
    else z += b(i + g() - 2) * x;
  }
  return blank ? BLANKf : z;    // 計算できたらzを返す
}

class mlogclass                 // mlogメインクラス
{
private:
  string fname;                 // ファイル名
  idd id;                       // データ識別情報
  ldata datalst;                // データリスト
  void put_resid(ostream& fout);
  // 解析結果のタイトルとデータ識別情報を出力する
  void put_stat(ostream& fout, const varclass& var); // 群ごとの統計量を出力する
  void put_mlogreg(ostream& fout, option& opt, const varclass& var);
  // 順序ロジスティック回帰分析結果を出力する
  void put_brci(ostream& fout, option& opt, const varclass& var) const;
  // 偏回帰係数と相対リスクの信頼区間を出力する
  void put_mlogtab(ostream& fout, const varclass& var);
  // 順序ロジスティック回帰分析の分散分析表を出力する
  int cal_mloglf(double& lf, const varclass& var);
  // 飽和モデルの尤度を求める
  void put_dns(option& opt, const varclass& var);
  // データファイル形式のロジットスコアと確率を出力する
  void put_dnsid(ostream& fout, const varclass& var);
  // DANSデータのタイトルとデータ識別情報を出力する
  void plot_logngp(option& opt, const varclass& var);
  // Ngraph用データ出力
  void plot_ngpid(ofstream& fout);
  // Ngraph用データのタイトルとデータ識別情報を出力する
  double cal_zscore(int bidx, int g, const vec<double>& bv, datac& cdata);
  // ロジットスコアを計算して返す
  void plot_score(ofstream& fout, const varclass& var);
  // Ngraph用ロジットスコアを出力する
  void plot_log(ofstream& fout, option& opt, const varclass& var);
  // Ngraph用ロジスティック曲線を出力する
public:
  mlogclass() : fname(), id(), datalst() // コンストラクタ
  {}
  ~mlogclass()                  // デストラクタ
  {}
  int get(const char* c_fname, option& opt); // データを読み込み累積する
  void put(option& opt);                     // 集計結果を出力する
};

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

void mlogclass::put_stat(ostream& fout, const varclass& var)
// 群ごとの統計量を出力する
{
  size_t dim = id.ni() + 1;       // 次数
  vecdbl meveci(dim, 0.0), vrveci(dim, 0.0);
  int i = 1;
  for (auto& gstati : var.gsmap())
  {
    int n = gstati.second.n();
    for (size_t j = 0; j < dim; ++j)
    {
      meveci[j] = gstati.second.smat(0, j);
      vrveci[j] = (n > 1) ? gstati.second.smat(j, j) / (double)(n - 1) : 0.0;
    }
    cout << "\n○順序" << i << ":" << id.item(0) << "=" << gstati.first;
    putmsdvec(cout, n, meveci, vrveci, id.ni() - 1);
    ++i;
  }
}

void mlogclass::put_mlogreg(ostream& fout, option& opt, const varclass& var)
// 順序ロジスティック回帰分析結果を出力する
{
  fout << "\n・最尤法による解析結果  反復回数=" << var.nrep();
  if (var.resvar() == reserr)
  { fout << ":Newton法発散！"; }
  else if (var.resvar() == resrmax)
  {
    if (var.nrep() == 1) fout << ":初期値を使用";
    else fout << ":最大反復回数";
  }
  fout << "\n  累積ロジットモデル:pk=1/{1+exp(-β0k-Σβj・xj)}\n";
  fout << "  pk:順序(k+1)以上の累積確率 β0k:kの切片 βj:変数xjの偏回帰係数\n";
  if (opt.res("v"))
  {
    fout << "  偏回帰係数初期値:";
    for (int i = 0; i < var.g() - 1; ++i)
    {
      if (i > 0)
      {
        if (i % 4 == 0) fout << "\n" << setw(19) << " ";
        else fout << " ";
      }
      fout << "β0" << i + 1 << "=" << var.b0(i);
    }
    for (int i = var.g() - 1; i < var.np(); ++i)
    {
      int j = i - (var.g() - 1);
      if (j % 4 == 0) fout << "\n" << setw(18) << " ";
      fout << " β" << i - var.g() + 2 << "=" << var.b0(i);
    }
    fout << "\n";
  }
  if (var.resvar() == reserr) return;
  fout << "                                            標準                   有意確率\n";
  fout << "変数  偏回帰係数   標準誤差    オッズ比  偏回帰係数  Waldのχ2乗       p値\n";
  constexpr int LENline = 78;   // 区切り線の長さ
  putline(fout, LENline);
  for (int i = 0; i < var.g() - 1; ++i)
  {
    double csqr = (var.vb(i) > 0.0) ? var.b(i) * var.b(i) / var.vb(i) : 0.0;
    ldbl pr = ctop(csqr, 1);
    fout << "切片" << setw(2) << i + 1 << setw(10) << var.b(i);
    fout << setw(11) << sqrt(var.vb(i)) << setw(36) << csqr;
    fout << setw(12) << pr << pind(pr) << "\n";
  }
  for (int i = var.g() - 1; i < var.np(); ++i)
  {
    double csqr = (var.vb(i) > 0.0) ? var.b(i) * var.b(i) / var.vb(i) : 0.0;
    double pr = ctop(csqr, 1);
    fout << "x" << setw(3) << i - var.g() + 2 << setw(12) << var.b(i);
    fout << setw(11) << sqrt(var.vb(i)) << setw(12) << exp(var.b(i));
    fout << setw(12) << var.bst(i);
    fout << setw(12) << csqr << setw(12) << pr << pind(pr) << "\n";
  }
  putline(fout, LENline);
  put_brci(fout, opt, var);
  if (opt.res("ltab")) put_mlogtab(fout, var);
}

void mlogclass::put_brci(ostream& fout, option& opt, const varclass& var) const
// 偏回帰係数と相対リスクの信頼区間を出力する
{
  fout << "\n変数  偏回帰係数" << setw(5) << opt.cc() * 100.0;
  fout << "%CI下限        上限    オッズ比" << setw(5) << opt.cc() * 100.0;
  fout << "%CI下限        上限\n";
  constexpr int LENline = 76;   // 区切り線の長さ
  putline(fout, LENline);
  double z = ptoz(1.0 - opt.cc());
  for (int i = 0; i < var.g() - 1; ++i)
  {
    double ci = z * sqrt(var.vb(i));
    double bl = var.b(i) - ci, bu = var.b(i) + ci;
    fout << "切片" << setw(2) << i + 1 << setw(10) << var.b(i);
    fout << setw(12) << bl << setw(12) << bu << "\n";
  }
  for (int i = var.g() - 1; i < var.np(); ++i)
  {
    double ci = z * sqrt(var.vb(i));
    double bl = var.b(i) - ci, bu = var.b(i) + ci;
    fout << "x" << setw(3) << i - var.g() + 2 << setw(12) << var.b(i);
    fout << setw(12) << bl << setw(12) << bu;
    fout << setw(12) << exp(var.b(i)) << setw(12) << exp(bl);
    fout << setw(12) << exp(bu) << "\n";
  }
  putline(fout, LENline);
}

void mlogclass::put_mlogtab(ostream& fout, const varclass& var)
// 順序ロジスティック回帰分析の分散分析表を出力する
{
  double l0 = 0.0;
  for (const auto& gstati : var.gsmap())
  {
    int nk = gstati.second.n();
    if (nk > 0 && nk < var.nt()) l0 += (double)nk * log((double)nk / (double)(var.nt()));
  }
  int dft = var.nt() - 1;
  double lf = 0.0;
  dft -= cal_mloglf(lf, var);
  int dfb = (int)(id.ni()) - 1;
  int dfl = dft - dfb;
  double r2 = 0.0;
  if (sgndiff(lf, l0) != 0) r2 = (var.lb() - l0) / (lf - l0);
  fout << "対数尤度:L(モデル)=" << var.lb() << "  L(切片)=" << l0 << "  飽和Lf=" << lf << "\n";
  fout << "擬似寄与率R2乗=" << r2 << "  AIC(赤池の情報量基準)=" << 2.0 * (dfb + 1 - var.lb()) << "\n";
  const string Head = "回帰とズレの検定";
  int w = strclen(Head, CHARwide);
  int npad = Head.size() - w;
  fout << "\n" << setw((68 + w) / 2 + npad) << Head;
  fout << "\n要因         -(対数尤度比)    自由度        χ2乗値   有意確率p値\n";
  constexpr int LENline = 68;
  putline(fout, LENline);
  double csqr = 2.0 * (var.lb() - l0);
  ldbl pr = ctop(csqr, dfb);
  fout << "回帰" << setw(22) << var.lb() - l0 << setw(10) << dfb;
  fout << setw(14) << csqr << setw(14) << pr << pind(pr) << "\n";
  csqr = 2.0 * (lf - var.lb());
  if (dfl < 1 || csqr < FDeps) csqr = 0.0;
  pr = ctop(csqr, dfl);
  fout << "ズレ(LOF)" << setw(17) << lf - var.lb() << setw(10) << dfl;
  fout << setw(14) << csqr << setw(14) << pr << pind(pr) << "\n";
  putline(fout, LENline);
  fout << "全体" << setw(22) << lf - l0 << setw(10) << dft << endl;
}

int mlogclass::cal_mloglf(double& lf, const varclass& var)
// 飽和モデルの尤度を求める
{
  vector<bool> dcvec(datalst.size(), true);
  int p = id.ni() - 1, na = var.g();
  vec<double> xiv(p), xjv(p);
  int nlf = 0, dci = 0;
  lf = 0.0;
  for (ldata_itr iitr = datalst.begin(); iitr != datalst.end(); ++iitr, ++dci)
  {
    if (dcvec[dci] == true && iitr->fdata(0) != BLANKf)
    {
      int ni = 1;
      vecint nivec(na, 0);
      int k = var.search_index(iitr->fdata(0));
      nivec[k] = 1;
      for (int j = 0; j < p; ++j) xiv.assign(j, iitr->fdata(j + 1));
      ldata_itr jitr = iitr;
      ++jitr;
      int dcj = dci + 1;
      while (jitr != datalst.end())
      {
        if (dcvec[dcj] == true && jitr->fdata(0) != BLANKf)
        {
          for (int j = 0; j < p; ++j) xjv.assign(j, jitr->fdata(j + 1));
          if (xiv == xjv)       // xが等しい
          {
            dcvec[dcj] = false; // 1度しか数えない
            ++ni;
            k = var.search_index(jitr->fdata(0));
            ++(nivec[k]);
          }
        }
        ++jitr;  ++dcj;
      }
      if (ni > 1)
      {
        nlf += ni - 1;
        for (int k = 0; k < na; ++k)
        {
          int nk = nivec[k];
          if (nk > 0 && nk < ni) lf += (double)nk * log((double)nk / (double)ni);
        }
      }
    }
  }
  return nlf;                   // 飽和モデルの重複数を返す(自由度減少数)
}

void mlogclass::put_dns(option& opt, const varclass& var)
// データファイル形式のロジットスコアと反応確率を出力する
{
  ofstream fout;
  if (ofopen(fout, (opt.dnsname()).c_str()) == false) return; // ファイルオープン失敗
  put_dnsid(fout, var);
  datac cdata(2 * var.g());
  for (auto& cdatai : datalst)
  {
    double y = cdatai.fdata(0);
    cdata.assign(0, y);
    if (y == BLANKf)
    { for (int i = 1; i < 2 * var.g(); ++i) cdata.assign(i, BLANKs); }
    else
    {
      for (int i = 1; i < var.g(); ++i) cdata.assign(i, var.zscore(i - 1, cdatai));
      double zk = cdata.fdata(1);
      if (zk == BLANKf)
      { for (int k = 0; k < var.g(); ++k) cdata.assign(var.g() + k, BLANKf); }
      else
      {
        double pk = 1.0 / (1.0 + exp(-zk));
        cdata.assign(var.g(), 1.0 - pk);
        double ppre = pk;
        for (int k = 2; k < var.g(); ++k)
        {
          zk = cdata.fdata(k);
          pk = 1.0 / (1.0 + exp(-zk));
          cdata.assign(var.g() + k - 1, ppre - pk);
          ppre = pk;
        }
        cdata.assign(2 * var.g() - 1, pk);
      }
    }
    cdata.put(fout);
  }
}

void mlogclass::put_dnsid(ostream& fout, const varclass& var)
// DANSデータのタイトルとデータ識別情報を出力する
{
  idd tid(id.title(), 2 * var.g());
  tid.assign_item(0, id.item(0));
  for (int i = 1; i < var.g(); ++i)
  {
    ostringstream ostr;
    ostr << "ロジットスコア" << i;
    tid.assign_item(i, ostr.str());
  }
  int i = 0;
  for (const auto& gstati : var.gsmap())
  {
    ostringstream ostr;
    ostr << id.item(0) << "=" << gstati.first << "になる確率";
    tid.assign_item(var.g() + i, ostr.str());
    ++i;
  }
  id.put_title(fout, formtag);
  id.put_stra(fout, formtag);
  tid.put_item(fout, formtag);
}

void mlogclass::plot_logngp(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 -w80 -c255,0,0 -zy-0.1,1.1,0.2\n";
    plot_ngpid(fout);
    plot_log(fout, opt, var);
    delml = true;
  }
  if (opt.ngp("score"))
  {
    if (delml) fout << NGPdelm << "\n";
    else fout << "-s1 -dmark,00 -o300 -zy-0.1,1.1,0.2\n";
    plot_ngpid(fout);
    plot_score(fout, var);
  }
}

void mlogclass::plot_ngpid(ofstream& fout)
// Ngraph用データのタイトルとデータ識別情報を出力する
{
  fout << "%";  puttitle(fout, Title);  fout << "%\n";
  id.put_title(fout, formngp);
  id.put_stra(fout, formngp);
  fout << "%目的変数y   :" << id.item(0) << "\n";
  for (size_t i = 1; i < id.ni(); ++i)
  { fout << "%説明変数x" << setw(3) << i << ":" << id.item(i) << "\n"; }
}

void mlogclass::plot_score(ofstream& fout, const varclass& var)
// Ngraph用ロジットスコアを出力する。
// プロット用に本来のロジットスコアの切片を切片1に固定する。
{
  for (int i = 1; i < var.g(); ++i)
  {
    if (i > 1) fout << NGPdelm << "\n";
    fout << "%ロジットスコア" << i << "\n";
    for (auto& cdatai : datalst)
    {
      double y = cdatai.fdata(0);
      if (y != BLANKf)
      {
        int k = var.search_index(y);
        y = (k < i) ? 0.0 : 1.0;
        double z = var.zscore(0, cdatai);
        if (z != BLANKf) fout << z << " " << y << "\n";
      }
    }
  }
}

void mlogclass::plot_log(ofstream& fout, option& opt, const varclass& var)
// Ngraph用ロジスティック曲線を出力する
// プロット用にロジットスコアの切片を切片1に固定し、x座標を定数の差分だけずらす。
{
  constexpr int RDdef = 50;     // ロジスティック曲線の分解数
  double z = 0.0, zmin = -5.0, zmax = 5.0, zmg = 0.1;
  for (auto& cdatai :  datalst)
  {
    double y = cdatai.fdata(0);
    if (y != BLANKf)
    {
      z = var.zscore(0, cdatai);
      if (z != BLANKf)
      {
        if (z < zmin) zmin = z;
        if (z > zmax) zmax = z;
      }
    }
  }
  double xdiff = var.b(0) - var.b(var.g() - 2);
  if (xdiff > 0.0) zmax += xdiff;
  else zmin += xdiff;
  if (opt.xmin() != BLANKf && opt.xmax() != BLANKf)
  { zmin = opt.xmin();  zmax = opt.xmax();  zmg = 0.0; }
  double range = zmax - zmin;
  double xstart = zmin - range * zmg, xend = zmax + range * zmg;
  double dx = (xend - xstart) / (double)RDdef;
  for (int i = 1; i < var.g(); ++i)
  {
    if (i > 1) fout << NGPdelm << "\n";
    fout << "%ロジスティック曲線" << i << "\n";
    xdiff = var.b(0) - var.b(i + var.g() - 4);
    double x = xstart;
    for (int i = 0; i <= RDdef; ++i, x += dx)
    {
      double pr = 1.0 / (1.0 + exp(-x));
      fout << x + xdiff << " " << pr << "\n";
    }
  }
}

int mlogclass::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 (id.ni() < 2) throw runtime_error("項目が2つ以上必要です！");
  return nc;
}

void mlogclass::put(option& opt) // 集計結果を出力する
{
  varclass var(id.ni());        // 変数オブジェクト
  var.cal(opt, datalst);        // 各種の計算実行
  if (opt.check(var.np()) != 0) return; // オプション指定エラー
  if (opt.res("res"))
  {
    put_resid(cout);
    if (opt.res("stat"))
    {
      put_stat(cout, var);
      cout << "\n○全体";
      putmsdvec(cout, var.nt(), var.mevec(), var.vrvec(), id.ni() - 1, 1);
      putmat(cout, "・相関行列(correlation coefficient matrix)", matr, var.nt(), var.smat(), id.ni() - 1, 1);
      if (var.resvar() == reserr)
      { cout << "\n※ロジスティック関数計算不可能！\n\n" << endl; }
    }
  }
  if (var.resvar() == reserr) return;
  var.mlog(opt, datalst);
  if (opt.res("res"))
  {
    if (opt.res("olra")) put_mlogreg(cout, opt, var);
    cout << "\n" << endl;
  }
  if (var.resvar() == reserr) return;
  if (opt.dnsname() != "") put_dns(opt, var);
  if (opt.ngpname() != "") plot_logngp(opt, var);
}

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