// disan.C Version 7.45  Copyleft (c) 杉本典夫
// Time-stamp: <2026-04-28 14:41:30 tac>
//
// 最初の項目を判別項目、以後の項目を説明変数とした逐次変数選択法による
// 判別分析を行う。
//
// 使用ライブラリ:libdans.so, libfun.so
// コンパイル法  :g++ -O -s -W disan disan.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 =            // プログラム使用法表示用
{
  "disan  Version 7.45          Copyleft(Vo\\oV) 2026 杉本典夫",
  "書式:disan [<オプション> …] [<データファイル> …]",
  "機能:最初の項目を判別項目、以後の項目を変数にした逐次変数選択法による",
  "     判別分析を行います。",
  "<オプション>",
  "  -cc[=<信頼係数>]:信頼係数を指定します。(0≦信頼係数<1)",
  "    -cc:-cc=0.95(オプション省略時)",
  "  -comp[={-|+}]:判別スコア計算時の欠損を補間します。(-comp:-comp=+)",
  "  -dns[={coef|score}]:DANSデータの種類を指定します。",
  "    coef:判別係数  score:判別スコア  -dns:-dns=score",
  "  -dnsname[=<出力ファイル名>]:DANSデータを出力するファイル名を指定します。",
  "    オプション省略時:DANSデータ出力せず  -dnsname:-dnsname=disan_dns.d",
  "  -f=<ファイル>:指定したファイルからオプションを読み込みます。",
  "  -fin=<取り込み基準>:変数の取り込み基準Finを指定します。(Fin≧0)",
  "  -fout=<追い出し基準>:変数の追い出し基準Foutを指定します。(0≦Fout≦Fin)",
  "  -g{1|2}[=<群コード>]:群1または群2のコードを指定します。",
  "  -h[=<コマンド名>]:オンラインマニュアルを表示します。(-h:-h=man)",
  "  -ix=<変数の開始番号>[-<終了番号>]x<変数の開始番号>[-<終了番号>][,…]",
  "    交互作用項目を作成する2つの変数の番号を指定します。",
  "  -ixc[={-|+}]:交互作用項目の変数を中心化するかどうか指定します。",
  "    -:中心化しない(オプション省略時)  +:中心化する  -ixc:ixc=+",
  "  -ngp={lci|log|score}[,…]:Ngraph用出力データを指定します。",
  "    lci:ロジスティック曲線の信頼区間  log:ロジスティック曲線",
  "    score:判別スコア  オプション省略時:全て出力",
  "  -ngpname[=<出力ファイル名>]:Ngraph用データの出力ファイル名を指定します。",
  "    オプション省略時:Ngraph用データ出力せず  -ngpname:-ngpname=disan_ngp.dat",
  "  -pi=<群1の事前確率>:群1の事前確率を指定します。(0＜事前確率＜1)",
  "  -pin=<取り込み基準>:変数の取り込み基準Pinを指定します。(0<Pin≦1)",
  "  -pout=<追い出し基準>:変数の追い出し基準Poutを指定します。(1≧Pout≧Pin)",
  "  -res=[-][{def|all|allx|selx|stata|statg}][,…]",
  "      :解析結果の出力項目を指定します。",
  "    -:-を付けた項目を出力しない(-res=-:全ての項目を出力しない)",
  "    def:オプション省略時と同じ項目を出力する  all:全ての項目を出力する",
  "    allx:全変数選択結果  selx:変数選択結果  stata:全体の基礎統計量",
  "    statg:群別基礎統計量",
  "  -x<変数番号>=<変数の値>:X軸にしない変数の値を指定します。",
  "  -xno=<変数の番号>:X軸にする変数の番号を指定します。",
  "  -xrange=<最小値>-<最小値>:X軸の最小値と最大値を指定します。"
};

const string Title = "=== 判別分析(discriminant analysis) ==="; // 表題
const string DNSNAMEdef = "disan_dns.d"; // 省略時DANSファイル名
const string NGPNAMEdef = "disan_ngp.dat"; // 省略時Ngraphファイル名
constexpr double FIOdef = 2.0;  // 省略時変数F値むの取り込み・追い出し基準
constexpr double PIOdef = 0.1;  // 省略時変数p値の取り込み・追い出し基準
constexpr int RDdef = 50;       // ロジスティック曲線の分解数省略時値

enum dnstype                    // DANSデータの種類
{
  dnsnon,                       // 出力なし
  dnscoef,                      // 判別スコア計算用判別係数
  dnsscore,                     // 判別スコア
  dnsetc                        // その他
};

class option                    // オプションクラス
{
private:
  string v_path;                // プログラムのパス
  string v_pname;               // プログラム名
  double v_cc;                  // 信頼係数
  bool v_comp;                  // 欠損補間の有無
  dnstype v_dns;                // DANS出力データの種類
  string v_dnsname;             // DANSファイル名
  pairdbl flevel;               // F値の取り込み基準と追い出し基準
  pairdbl plevel;               // P値の取り込み基準と追い出し基準
  string v_ngpname;             // Ngraphファイル名
  lstr ngplst;                  // Ngraph出力リスト
  mapsb ngpmap;                 // Ngraph出力の有無マップ
  vecstr gvec;                  // 群コード
  lstr ixstrlst;                // 交互作用項目の変数リスト
  lino v_ixlst;                 // 交互作用項目の変数ペアリスト
  bool v_ixc;                   // 交互作用項目の中心化フラッグ
  double v_pi;                  // ロジスティック曲線における群1の事前確率
  lstr reslst;                  // 解析結果出力リスト
  mapsb resmap;                 // 解析結果出力の有無マップ
  lstr xlst;                    // 変数リスト
  vecdbl xvec;                  // 変数ベクトル
  int v_xno;                    // X軸にする変数の番号
  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; }
  bool comp() const             // 欠損補間の有無を返す
  { return v_comp; }
  dnstype dns() const           // DANS出力データの種類を返す
  { return v_dns; }
  const string& dnsname() const // DANSデータファイル名を返す
  { return v_dnsname; }
  double fin() const            // F値の取り込み基準を返す
  { return flevel.first; }
  double fout() const           // F値の追い出し基準を返す
  { return flevel.second; }
  const string& g(int idx) const // 群コードを返す
  { return gvec[idx]; }
  const lino& ixlst()    // 交互作用項目の変数番号ペアリストを返す
  { return v_ixlst; }
  bool ixc()                   // 交互作用項目の中心化フラッグを返す
  { return v_ixc; }
  double pi() const   // ロジスティック曲線における群1の事前確率を返す
  { return v_pi; }
  double pin() const            // P値の取り込み基準を返す
  { return plevel.first; }
  double pout() const           // P値の追い出し基準を返す
  { return plevel.second; }
  const string& ngpname() const // Ngraphデータファイル名を返す
  { return v_ngpname; }
  bool ngp(const string& code); // Ngraph出力フラッグを返す
  bool res(const string& code); // 分析結果出力フラッグを返す
   double x(size_t xno) const    // 変数の値を返す
  { return (0 < xno && xno <= xvec.size()) ? xvec[xno - 1] : BLANKf; }
  int xno() const               // X軸にする変数の番号を返す
  { return v_xno; }
  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オプションのセット
  int set_ix(int nx);           // -ixオプションのセット
  int set_x(int nx);            // -xオプションのセット
};

option::option(const char* c_pname) : v_path(), v_pname(), v_cc(CCdef), v_comp(false), v_dns(dnsscore), v_dnsname(), flevel(FIOdef, FIOdef), plevel(1.0, 1.0), gvec(2), ixstrlst(), v_ixlst(), v_ixc(false), v_ngpname(), v_pi(0.0), reslst(), resmap(), xlst(), xvec(), v_xno(0), xrange(BLANKf, BLANKf)
// コンストラクタ
{
  if (c_pname == 0) return;     // 初期設定せず
  v_pname = getname(c_pname);   // パス名なしのプログラム名取得
  v_path = c_pname;
  strsize_t idx = v_path.find(v_pname);
  if (idx != string::npos) v_path.erase(idx); // プログラム名削除
  gvec[0] = "1";
}

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] == "comp")    // 欠損補完の有無の指定
  {
    if (opt[1].empty() || opt[1] == "+") v_comp = true;
    else if (opt[1] == "-") v_comp = false;
    else err = badopt(pname(), optstr, USAGE);
  }
  else if (opt[0] == "dns")     // DANSデータの種類指定
  {
    if (opt[1].empty() || opt[1] == "score") v_dns = dnsscore;
    else if (opt[1] == "coef") v_dns = dnscoef;
    else err = badopt(pname(), optstr, USAGE);
  }
  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] == "f")       // オプション指定ファイル
  { err = fgetopt(*this, opt[1], true); }
  else if (opt[0] == "fin")     // F値の取り込み基準指定
  {
    flevel.first = opt[1].empty() ? FIOdef : strtodf(opt[1]);
    if (flevel.first < 0.0) err = badopt(pname(), optstr, USAGE, 0.0, FDmax);
    if (flevel.second > flevel.first) flevel.second = flevel.first;
  }
  else if (opt[0] == "fout")    // F値の追い出し基準指定
  {
    flevel.second = opt[1].empty() ? flevel.first : strtodf(opt[1]);
    if (flevel.second < 0.0 || flevel.first < flevel.second) err = badopt(pname(), optstr, USAGE, 0.0, flevel.first);
  }
  else if (opt[0] == "g1")      // 群1のコード指定
  { gvec[0] = opt[1]; }
  else if (opt[0] == "g2")      // 群2のコード指定
  { gvec[1] = opt[1]; }
  else if (opt[0] == "h")       // オンラインヘルプ表示
  { err = help(pname(), USAGE, opt[1]); }
  else if (opt[0] == "ix")      // 交互作用項目の変数ペア指定
  {
    if (strsplit(ixstrlst, opt[1], CHRdelm) < 1) err = badopt(pname(), optstr, USAGE);
  }
  else if (opt[0] == "ixc")     // 交互作用項目の変数の中心化指定
  {
    if (opt[1].empty() || opt[1] == "+" ) v_ixc = true;
    else if (opt[1] == "-") v_ixc = false;
    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] == "pi")      // ロジスティック曲線における群1の事前確率指定
  {
    v_pi = opt[1].empty() ? 0.0 : strtodf(opt[1]);
    if (v_pi < 0.0 || 1.0 <= v_pi) err = badopt(pname(), optstr, USAGE, 0.0, 1.0);
  }
  else if (opt[0] == "pin")     // p値の取り込み基準指定
  {
    plevel.first = opt[1].empty() ? PIOdef : strtodf(opt[1]);
    if (plevel.first <= 0.0 || 1.0 < plevel.first) err = badopt(pname(), optstr, USAGE, 0.0, 1.0);
  }
  else if (opt[0] == "pout")    // p値の追い出し基準指定
  {
    plevel.second = opt[1].empty() ? plevel.first : strtodf(opt[1]);
    if (plevel.second < plevel.first || 1.0 < plevel.second) err = badopt(pname(), optstr, USAGE, plevel.first, 1.0);
  }
  else if (opt[0] == "res")     // 分析結果の出力指定
  {
    if (strsplit(reslst, opt[1], CHRdelm) < 1) err = badopt(pname(), optstr, USAGE);
  }
  else if (opt[0] == "xno")     // X軸にする変数の番号指定
  {
    char* ptr = nullptr;
    v_xno = strtol(opt[1].c_str(), &ptr, 0);
    if (v_xno < 1) err = badopt(pname(), optstr, USAGE, 1.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 if (opt[0][0] == 'x')    // 変数の値指定
  {
    if (strsplit(xlst, optstr, CHRdelm) < 1) err = badopt(pname(), optstr, USAGE);
  }
  else                          // 間違ったオプション指定
  { err = nonopt(pname(), optstr, USAGE); }
  if (err > 0) throw help_except("ヘルプ表示");
  if (err < 0) throw option_error(string("オプション指定違反→") + optstr);
}

int option::check(size_t ni)    // オプションのセットとチェック
{
  if (ni < 2) return disperr(pname(), "項目数が足りません！");
  if ((gvec[0]).empty() && (gvec[1]).empty()) return disperr(pname(), "群コードが指定されていません！");
  int nx = (int)ni - 1;         // 交互作用項目を含まない変数の数
  if (nx < v_xno) return disperr(pname(), "X軸にする変数の番号が間違っています！");
  if (set_ngp() < 0) return -2; // -ngpオプションのエラー
  if (set_res() < 0) return -3; // -resオプションのエラー
  int nxix = set_ix(nx);        // -ixオプションのセット
  if (nxix < 0) return -4;      // 交互作用項目を含んだ変数の数
  if (set_x(nx) < 0) return -5; // -xオプションのセット
  return 0;                     // エラーなし
}

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

int option::set_ix(int nx)      // -ixオプションのセット
{
  if (ixstrlst.empty()) return nx;
  int ixno = nx + 1;            // 交互作用項目の最初の番号
  for (const string& ixstr : ixstrlst)
  {
    string opt[2];
    if (getopt(opt, ixstr, 'x') > 0)
    {
      int low0 = 0, up0 = 0, low1 = 0, up1 = 0;
      if (getrange(low0, up0, opt[0]) < 1)
      { return disperr(pname(), "1番目の変数番号が間違っています→", ixstr); }
      else
      {
        if (getrange(low1, up1, opt[1]) < 1)
        { return disperr(pname(), "2番目の変数番号が間違っています→", ixstr); }
        else
        {
          if (low0 < 1 || low0 > up0) return disperr(pname(), "1番目の変数番号が間違っています→", ixstr);
          if (low1 < 1 || low1 > up1) return disperr(pname(), "2番目の変数番号が間違っています→", ixstr);
          if (ixno <= up0 || ixno <= up1) return disperr(pname(), "変数番号が間違っています→", ixstr);
          for (int j = low1; j <= up1; ++j) // 2番目の番号
          {
            for (int i = low0; i <= up0; ++i) // 1番目の番号
            {
              pairint ixpi {i, j};
              v_ixlst.push_back(ixpi);
            }
          }
          ixno += (up0 - low0 + 1) * (up1 - low1 + 1);
        }
      }
    }
  }
  return ixno - 1;              // 交互作用項目を含めた変数の数
}

int option::set_x(int nx)       // -xオプションのセット
{
  if (xlst.empty()) return 0;   // 変数指定なし
  if ((int)(xlst.size()) > nx) return disperr(pname(), "変数の値が多すぎます！");
  xvec.resize(nx);
  xvec.assign(nx, BLANKf);
  for (const string& xstr : xlst)
  {
    string opt[2];
    if (getopt(opt, xstr, OPTdev) > 0)
    {
      char* ptr = nullptr;
      int xno = strtol(opt[0].substr(1).c_str(), &ptr, 0);
      if (xno < 1 || nx < xno || xno == v_xno) return disperr(pname(), "変数の項目番号が間違っています→", xstr);
      double fdata = strtodf(opt[1]);
      if (fdata == BLANKf) return disperr(pname(), "変数の値が間違っています→", xstr);
      xvec[xno - 1] = fdata;
    }
  }
  return (int)(xlst.size());    // 変数の数を返す
}

class varclass                 // 変数クラス
{
private:
  vector<mat<double>> smatvec;  // 群別積和行列ベクトル
  vecint nvec;                  // 群別例数ベクトル
  mat<double> v_dmat;           // 群内積和行列
  vecdbl v_mevec;               // 総平均値ベクトル
  vecdbl midvec;                // 2群の平均値の中央値ベクトル
  vecdbl v_vrvec;               // 総群内分散ベクトル
  vecint checkvec;              // 変数選択用チェックベクトル
  int v_nx;                     // 変数の数
  int v_nxix;                   // 交互作用項目を含んだ変数の数
  int v_nt;                     // 総例数
  restype v_resvar;             // 処理結果の種類
  int set_mean(vecdbl& mevec, option& opt, ldata& datalst);
  // 交互作用項計算用平均値ベクトルをセットする
  int set_fdata(vecdbl& fdatavec, option& opt, const vecdbl& mevec, datac& cdata) const;
  // 交互作用項目を含むデータをセットする
  void set_dmat();              // 変数選択用積和行列を準備する
public:
  varclass(size_t ni, option& opt); // コンストラクタ
  ~varclass()                   // デストラクタ
  {}
  int n(int k) const            // 群別例数を返す
  { return (0 <= k && k < 2) ? nvec[k] : 0; }
  double me(int k, int j) const // 群別平均値を返す
  { return (0 <= k && k < 2 && 0 <= j && j <= np()) ? (smatvec[k])(0, j) : 0.0; }
  double vr(int k, int j) const // 群別分散を返す
  { return (0 <= k && k < 2 && 0 <= j && j <= np()) ? (smatvec[k])(j, j) : 0.0; }
  double mid(int idx) const     // 2群の平均値の中央値を返す
  { return (0 <= idx && idx < np()) ? midvec[idx] : 0.0; }
  double vrt(int idx) const     // 総群内分散を返す
  { return (0 <= idx && idx < np()) ? v_vrvec[idx] : 0.0; }
  const mat<double>& dmat() const // 群内積和行列を返す
  { return v_dmat; }
  const vecdbl& mevec() const   // 全体平均値ベクトルを返す
  { return v_mevec; }
  const vecdbl& vrvec() const   // 全体分散ベクトルを返す
  { return v_vrvec; }
  int nx() const                // 変数の数を返す
  { return v_nx; }
  int nxix() const              // 交互作用項目を含んだ変数の数を返す
  { return v_nxix; }
  int np() const                // パラメータ数を返す
  { return nxix() + 1; }
  int nt() const                // 総例数を返す
  { return v_nt; }
  int df1() const               // 第1自由度を返す
  { return checkvec[0]; }
  int df2() const               // 第2自由度を返す
  { return nt() - 1 - df1(); }
  int dft() const               // 総自由度を返す
  { return nt() - 2; }
  int check(int idx) const      // 変数選択フラッグを返す
  { return (0 <= idx && idx < np()) ? checkvec[idx] : 0; }
  double a(int idx) const         // 判別係数を返す
  { return (0 <= idx && idx <= np()) ? v_dmat(idx, np()) : 0.0; }
  double va(int idx) const        // 判別係数の分散を返す
  { return (0 <= idx && idx <= np()) ? v_dmat(idx, idx) : 0.0; }
  double d2() const             // 判別残差を返す
  { return -v_dmat(np(), np()); }
  restype resvar() const        // 処理結果の種類を返す
  { return v_resvar; }
  void cal(option& opt, ldata& datalst); // データを集計して各種統計量を求める
  void stepsel(double fin, double fout, double pin, double pout);
  // 変数増減法による判別関数の計算
  int sel_var(bool in, double fin, double fout, double pin, double pout, double w);
  // 変数増減法による変数の選択
  int calscore(double& z, option& opt, datac& cdata) const;
  // 判別スコアを計算し、群インデックスを返す
  double calrange(double& xmin, double& xmax, option& opt, datac& cdata) const;
  // ロジスティック曲線の最小値と最大値を計算する
};

varclass::varclass(size_t ni, option& opt) : smatvec(), nvec(), v_dmat(), v_mevec(), midvec(), v_vrvec(), checkvec(), v_nx(0), v_nxix(0), v_nt(0), v_resvar(resnon)
// コンストラクタ
{
  v_nx = ni - 1;
  v_nxix = v_nx + (opt.ixlst()).size();
  size_t dim = v_nxix + 2;
  smatvec.resize(2);  nvec.resize(2);
  (smatvec[0]).resize(dim, dim);  (smatvec[1]).resize(dim, dim);
  v_dmat.resize(dim, dim);
  smatvec[0] = 0.0;  smatvec[1] = 0.0;  nvec.assign(2, 0);
  v_mevec.resize(dim - 1, 0.0);  midvec.resize(dim - 1, 0.0);
  v_vrvec.resize(dim - 1, 0.0);  checkvec.resize(dim - 1, 0);
}

void varclass::cal(option& opt, ldata& datalst)
// データを集計して各種統計量を求める
{
  vecdbl mevec(nx() + 1);       // 交互作用項計算用平均値ベクトル
  v_nt = set_mean(mevec, opt, datalst);
  if (v_nt < np())              // 判別分析不可能
  {
    v_resvar = reserr;
    return;
  }
  int dim = np() + 1;
  vecdbl fdatavec(dim, 0.0);
  fdatavec[0] = 1.0;
  for (auto& cdatai : datalst)
  {
    int gidx = set_fdata(fdatavec, opt, mevec, cdatai);
    if (gidx == 0 || gidx == 1) // 群1と2の欠損無ケースだけ集計する
    {
      for (int i = 0; i < np(); ++i)
      {
        double x = fdatavec[i];
        for (int j = i; j < dim; ++j) smatvec[gidx].add(i, j, x * fdatavec[j]);
      }
    }
  }
  v_dmat = smatvec[0] + smatvec[1]; // 全体の単純積和行列
  v_dmat.sweep(0);
  for (int i = 0; i < np(); ++i) v_mevec[i] = v_dmat(0, i); // 総平均値
  nvec[0] = (int)(smatvec[0](0, 0));  nvec[1] = (int)(smatvec[1](0, 0));
  if (nvec[0] < 1 || nvec[1] < 1) // 判別分析不可能
  {
    v_resvar = reserr;
    return;
  }
  for (int k = 0; k < 2; ++k)
  {
    if (nvec[k] > 0)
    {
      for (int i = 0; i < np(); ++i)
      {
        for (int j = i + 1; j < dim; ++j) smatvec[k].assign(j, i, smatvec[k](i, j));
      }
      smatvec[k].sweep(0);
    }
  }
  v_dmat = smatvec[0] + smatvec[1]; // 群内積和行列
  set_dmat();                       // 変数選択用積和行列を準備する
  v_resvar = resok;
}

int varclass::set_mean(vecdbl& mevec, option& opt, ldata& datalst)
// 交互作用項計算用平均値ベクトルをセットする
{
  int ni = nx() + 1;            // 項目数
  mevec.assign(ni, 0.0);
  vecdbl fdatavec(ni, 0.0);
  int nt = 0;
  for (auto& cdatai : datalst)
  {
    if (cdatai.sdata(0) != BLANKs)
    {
      int gidx = -1;
      if (opt.g(0).empty())     // 群1コード指定なし
      { gidx = (cdatai.sdata(0) == opt.g(1)) ? 1 : 0; }
      else if (cdatai.sdata(0) == opt.g(0))
      { gidx = 0; }
      else if (opt.g(1).empty() || cdatai.sdata(0) == opt.g(1))
      { gidx = 1; }
      if (gidx < 0)             // 群1と2以外は集計から除外する
      { cdatai.assign(0, BLANKs); } // 後の集計のために群コードを欠損にする
      else
      {
        bool blank = false;
        int j = 1;
        while (blank == false && j < ni)
        {
          fdatavec[j] = cdatai.fdata(j);
          if (fdatavec[j] == BLANKf) blank = true;
          else ++j;
        }
        if (blank == false)     // 完全揃いケースだけ集計する
        {
          ++nt;  fdatavec[0] = gidx;
          for (int i = 0; i < ni; ++i) mevec[i] += fdatavec[i];
        }
      }
    }
  }
  for (int j = 0; j < ni; ++j) mevec[j] /= (double)nt;
  return nt;                    // 総例数を返す
}

int varclass::set_fdata(vecdbl& fdatavec, option& opt, const vecdbl& mevec, datac& cdata) const
// 交互作用項目を含むデータをセットする
{
  if (cdata.sdata(0) == BLANKs) return -1; // 群コード欠損
  int gidx = -1;
  if (opt.g(0).empty())         // 群1コード指定なし
  { gidx = (cdata.sdata(0) == opt.g(1)) ? 1 : 0; }
  else if (cdata.sdata(0) == opt.g(0))
  { gidx = 0; }
  else if (opt.g(1).empty() || cdata.sdata(0) == opt.g(1))
  { gidx = 1; }
  int ni = nx() + 1;            // 項目数
  bool blank = false;
  for (int j = 1; j < ni; ++j)
  {
    fdatavec[j] = cdata.fdata(j);
    if (fdatavec[j] == BLANKf) blank = true;
  }
  if ((opt.ixlst()).empty()) return blank ? gidx + 2 : gidx; // 交互作用項目なし
  int j = ni;
  for (const auto& linoi : opt.ixlst())
  {
    double x1 = fdatavec[linoi.first], x2 = fdatavec[linoi.second];
    if (x1 == BLANKf || x2 == BLANKf)
    { fdatavec[j] = BLANKf;  blank = true; }
    else
    {
      if (opt.ixc())            // 偏差に変換する
      {
        if (linoi.first < (int)(cdata.size())) x1 -= mevec[linoi.first];
        if (linoi.second < (int)(cdata.size())) x2 -= mevec[linoi.second];
      }
      fdatavec[j] = x1 * x2;
    }
    ++j;
  }
  return blank ? gidx + 2 : gidx; // 群インデックスを返す(0/1:欠損無 2/3:欠損有)
}

void varclass::set_dmat()       // 変数選択用積和行列を準備する
{
  midvec[0] = (double)(nt());
  v_dmat *= 1.0 / (double)(dft()); // 群内分散に変換する
  v_dmat.assign(0, 0, 0.0);  v_dmat.assign(0, np(), 0.0);
  v_dmat.assign(np(), 0, 0.0);  v_dmat.assign(np(), np(), 0.0);
  v_vrvec[0] = v_dmat(0, 0);
  for (int i = 1; i < np(); ++i)
  {
    midvec[i] = (smatvec[0](0, i) + smatvec[1](0, i)) * 0.5;
    v_dmat.assign(0, i, ((smatvec[0])(0, i) + (smatvec[1])(0, i)) * 0.5);
    v_dmat.assign(i, 0, -v_dmat(i, 0));
    v_dmat.assign(i, np(), (smatvec[0])(0, i) - (smatvec[1])(0, i));
    v_dmat.assign(np(), i, v_dmat(i, np()));
    v_vrvec[i] = v_dmat(i, i);  // 群内分散
  }
}

void varclass::stepsel(double fin, double fout, double pin, double pout)
// 変数増減法による判別関数の計算
{
  if (fin <= 0.0 && pin >= 1.0) // 全変数選択
  {
    for (int i = 1; i <= nxix(); ++i)
    {
      if (checkvec[i] == 0) { v_dmat.sweep(i);  checkvec[i] = 1; }
    }
    checkvec[0] = nxix();
  }
  else
  {
    int i = 0;
    bool in = true;
    double w = (double)(dft()) / (double)(nvec[0]) * (double)(nt()) / (double)(nvec[1]);
    do {
      in = true;
      while ((i = sel_var(in, fin, fout, pin, pout, w)) > 0)
      {
        v_dmat.sweep(i);  checkvec[i] = in ? 1 : 0;
        checkvec[0] += in ? 1 : -1;
        in = !in;
      }
    } while (in == false);        // 取り込む変数がなくなるまで続ける
  }
}

int varclass::sel_var(bool in, double fin, double fout, double pin, double pout, double w)
// 変数増減法による変数の選択
{
  int idx = 0, inout = in ? 1 : 0;
  double dsel = in ? FDmin : FDmax;
  bool fsel = true;             // F値による変数選択
  if (pin < 1.0) fsel = false;  // 有意確率pによる変数選択
  for (int j = 1; j < np(); ++j)
  {
    if (checkvec[j] != inout && fabs(v_dmat(j, j)) > FDeps)
    {
      double dj = v_dmat(j, np()) * v_dmat(j, np()) / v_dmat(j, j);
      if ((in && dj > dsel) || (in == false && dj < dsel))
      { dsel = dj;  idx = j; }
    }
  }
  if (idx < 1) return idx;      // 選択すべき変数なし
  if (in)                       // 変数取り込み
  {
    int df = nt() - checkvec[0] - 2;
    double f = dsel * (double)df / (w - v_dmat(np(), np()));
    if (fsel)
    {
      if (f <= fin) idx = 0;
    }
    else
    {
      ldbl pr = ftop(f, 1, df);
      if (pr >= (ldbl)pin) idx = 0;
    }
  }
  else                          // 変数追い出し
  {
    int df = nt() - checkvec[0] - 1;
    double f = dsel * (double)df / (w - v_dmat(np(), np()) - dsel);
    if (fsel)
    {
      if (f > fout) idx = 0;
    }
    else
    {
      ldbl pr = ftop(f, 1, df);
      if (pr < (ldbl)pout) idx = 0;
    }
  }
  return idx;                   // 選択すべき変数のインデックスを返す
}

int varclass::calscore(double& z, option& opt, datac& cdata) const
// 判別スコアを計算し、群インデックスを返す
{
  vecdbl fdatavec(np(), 0.0);
  int gidx = set_fdata(fdatavec, opt, v_mevec, cdata);
  if (gidx < 0) return gidx;            // 群コード不適切で計算できず
  if (gidx > 1 && opt.comp() == false)  // 欠損ありで補完なし
  {
    z = BLANKf;
    return -1;
  }
  z = v_dmat(0, np());
  bool blank = false;
  for (int i = 1; blank == false && i < np(); ++i)
  {
    if (checkvec[i])
    {
      double xi = fdatavec[i];
      if (xi == BLANKf) xi = midvec[i]; // 2群の平均値の中央値で補完する
      z += v_dmat(i, np()) * xi;
    }
  }
  return gidx;                  // 群インデックスを返す
}

double varclass::calrange(double& xmin, double& xmax, option& opt, datac& cdata) const
// ロジスティック曲線の最小値と最大値を計算する
{
  if (opt.xmin() != BLANKf && opt.xmax() != BLANKf) // X軸の範囲指定あり
  {
    xmin = opt.xmin();  xmax = opt.xmax();
    return xmax - xmin;         // 範囲を返す
  }
  double pr = 0.001;            // 最小確率
  double zmin = log(pr / (1 - pr));
  pr = 1 - pr;                  // 最大確率
  double zmax = log(pr / (1 - pr));
  if (opt.xno() > 0 && checkvec[opt.xno()]) // X軸にする変数が指定されている時
  {
    cdata.assign(opt.xno(), 0.0);
    double z = 0.0;
    calscore(z, opt, cdata);
    double axno = v_dmat(opt.xno(), np());
    xmin = (zmin - z) / axno;
    xmax = (zmax - z) / axno;
    if (axno < 0.0)             // 判別係数が負の時
    {
      z = xmin;
      xmin = xmax;
      xmax = z;
    }
  }
  else
  { xmin = zmin;  xmax = zmax; }
  return xmax - xmin;           // 範囲を返す
}

class disanclass                // disanメインクラス
{
private:
  string fname;                 // ファイル名
  idd id;                       // データ識別情報
  ldata datalst;                // データリスト
  void put_resid(ostream& fout, option& opt);
  // 解析結果のタイトルとデータ識別情報を出力する
  int set_itemix(idd& tid, option& opt, int ni) const;
  // 交互作用項目名をセットする
  void put_statg(ostream& fout, option& opt, const varclass& var);
  // 群別統計量を出力する
  void put_distab(ostream& fout, option& opt, double fin, double flout, double pin, double pout, const varclass& var);
  // 判別分析結果を出力する
  void put_arci(ostream& fout, option& opt, const varclass& var) const;
  // 判別係数の信頼区間を出力する
  void put_hotest(ostream& fout, const varclass& var) const;
  // HotellingのT2乗検定結果と推定結果を出力する
  void put_dnsid(ofstream& fout, option& opt);
  // DANSデータのタイトルとデータ識別情報を出力する
  void put_coef(ofstream& fout, const varclass& var);
  // 判別スコア計算用判別係数をデータファイル形式で出力する
  void put_score(ofstream& fout, option& opt, const varclass& var);
  // 判別スコアと群1に属す確率をデータファイル形式で出力する
  void plot_disngp(option& opt, const varclass& var);
  // Ngraph用データ出力
  void plot_ngpid(ofstream& fout, option& opt, const varclass& var);
  // Ngraph用データのタイトルとデータ識別情報を出力する
  void plot_log(ofstream& fout, option& opt, const varclass& var);
  // Ngraph用ロジスティック曲線を出力する
  void plot_lci(ofstream& fout, option& opt, const varclass& var);
  // Ngraph用ロジスティック曲線の信頼区間を出力する
  void plot_score(ofstream& fout, option& opt, const varclass& var);
  // Ngraph用判別スコアを出力する
public:
  disanclass() : fname(), id(), datalst() // コンストラクタ
  {}
  ~disanclass()                 // デストラクタ
  {}
  int get(const char* c_fname, option& opt); // データを読み込み累積する
  void put(option& opt);                     // 集計結果を出力する
};

void disanclass::put_resid(ostream& fout, option& opt)
// 解析結果のタイトルとデータ識別情報を出力する
{
  puttitle(fout, Title);
  fout << "\n";
  id.put_title(fout);
  id.put_stra(fout);
  idd tid(id);
  tid.set(id.title(), id.ni() + (opt.ixlst()).size()); // 項目数を変更する
  set_itemix(tid, opt, id.ni());
  fout << "判別項目:" << tid.item(0) << "\n";
  for (size_t i = 1; i < tid.ni(); ++i)
  {
    fout << "変数x" << setw(3) << i;
    fout << ":" << tid.item(i) << "\n";
  }
}

int disanclass::set_itemix(idd& tid, option& opt, int ni) const
// 交互作用項目名をセットする
{
  ostringstream ostr;
  int j = ni;
  for (const auto& linoi : opt.ixlst())
  {
    ostr.str("");
    if (opt.ixc())              // 偏差変換
    {
      if (linoi.first < ni) ostr << tid.item(linoi.first) << "偏差 x ";
      else ostr << tid.item(linoi.first) << " x ";
      if (linoi.second < ni) ostr << tid.item(linoi.second) << "偏差";
      else ostr << tid.item(linoi.second);
    }
    else
    { ostr << tid.item(linoi.first) << " x " << tid.item(linoi.second); }
    tid.assign_item(j, ostr.str());
    ++j;
  }
  return ni - j;                // セットした交互作用項目名数を返す
}

void disanclass::put_statg(ostream& fout, option& opt, const varclass& var)
// 群別統計量を出力する
{
  vecdbl mkvec(var.np(), 0.0), vkvec(var.np(), 0.0);
  for (int k = 0; k < 2; ++k)
  {
    if (var.n(k) > 0)
    {
      for (int j = 0; j < var.np(); ++j)
      {
        mkvec[j] = var.me(k, j);
        vkvec[j] = (var.n(k) > 1) ? var.vr(k, j) / (double)(var.n(k) - 1) : 0.0;
      }
      fout << "\n○群" << k + 1 << ":" << id.item(0) << "=";
      if (opt.g(k).empty()) fout << opt.g((k == 0) ? 1 : 0) << " 以外";
      else fout << opt.g(k);
      putmsdvec(fout, var.n(k), mkvec, vkvec);
    }
  }
}

void disanclass::put_distab(ostream& fout, option& opt, double fin, double flout, double pin, double pout, const varclass& var)
// 判別分析結果を出力する
{
  if (fin <= 0.0 && pin >= 1.0)
  { fout << "\n・全変数を選択した結果\n"; }
  else
  {
    fout << "\n・前進的変数増減法(stepwise forward selection method)による変数選択結果\n";
    if (pin < 1.0) fout << "  取り込み基準:p値＜Pin=" << pin << "  追い出し基準:p値≧Pout=" << pout << "\n";
    else fout << "  取り込み基準:偏F値＞Fin=" << fin << "  追い出し基準:偏F値≦Fout=" << flout << "\n";
  }
  fout << "                                             標準                  有意確率\n";
  fout << "変数    判別係数   標準誤差 exp(判別係数)  判別係数       偏F値       p値\n";
  constexpr int LENline = 78;   // 区切り線の長さ
  putline(fout, LENline);
  double w = (double)(var.n(0)) * (double)(var.n(1)) / (double)(var.nt());
  fout << "切片" << setw(12) << var.a(0) << "\n";
  for (int i = 1; i < var.np(); ++i)
  {
    if (var.check(i))
    {
      double ai = var.a(i);
      double di = (sgn(var.va(i)) == 0) ? 0.0 : ai * ai / var.va(i);
      double f = (var.df2() > 0) ? (double)(var.df2()) * di / ((double)(var.dft()) / w + var.d2() - di) : 0.0;
      double va = (sgn(f) == 0) ? 0.0 : ai * ai / f;
      ldbl pr = ftop(f, 1, var.df2());
      fout << "x" << setw(3) << i << setw(12) << ai;
      fout << setw(11) << sqrt(va) << setw(12) << exp(ai);
      fout << setw(12) << ai * sqrt(var.vrt(i)) << setw(12) << f;
      fout << setw(12) << pr << pind(pr) << "\n";
    }
  }
  putline(fout, LENline);
  if (var.df1() < 1) return;    // 選択変数なし
  put_arci(fout, opt, var);
  put_hotest(fout, var);
}

void disanclass::put_arci(ostream& fout, option& opt, const varclass& var) const
// 判別係数の信頼区間を出力する
{
  fout << "\n変数    判別係数" << setw(5) << opt.cc() * 100.0;
  fout << "%CI下限        上限 exp(判別係数)" << setw(3) << opt.cc() * 100.0;
  fout << "%CI下限        上限\n";
  constexpr int LENline = 76;   // 区切り線の長さ
  putline(fout, LENline);
  double w = (double)(var.n(0)) * (double)(var.n(1)) / (double)(var.nt());
  for (int i = 1; i < var.np(); ++i)
  {
    if (var.check(i))
    {
      double ai = var.a(i);
      double di = (sgn(var.va(i)) == 0) ? 0.0 : ai * ai / var.va(i);
      double va = (var.df2() > 0) ? var.va(i) * ((double)(var.dft()) / w + var.d2() - di) / (double)(var.df2()) : 0.0;
      double ci = (var.df2() > 0) ?  calmcci(mcfis, opt.cc(), 2, 1, var.df2(), 1, va) : 0.0;
      if (ci > FDmax) ci = 0.0;
      double al = ai - ci, au = ai + ci;
      fout << "x" << setw(3) << i << setw(12) << ai;
      fout << setw(12) << al << setw(12) << au;
      fout << setw(12) << exp(ai) << setw(12) << exp(al);
      fout << setw(12) << exp(au) << "\n";
    }
  }
  putline(fout, LENline);
}

void disanclass::put_hotest(ostream& fout, const varclass& var) const
// HotellingのT2乗検定結果と推定結果を出力する
{
  double w = (double)(var.n(0)) * (double)(var.n(1)) / (double)(var.nt());
  double cf = (var.nt() > 2) ? (double)(var.df2()) / (double)(var.dft()) / (double)(var.df1()) : 0.0;
  double t2 = w * var.d2();
  double f = t2 * cf;
  double r2 = f * (double)(var.df1()) / (f * (double)(var.df1()) + (double)(var.df2()));
  double z = sqrt(var.d2()) * 0.5;
  ldbl pr = sndis(z);
  fout << "・判別効率(Mahalanobisの平方距離):D2乗=" << var.d2() << "  相関比:η=" << sqrt(r2) << "\n";
  fout << "  D/2(正規分布z値)=" << z << "  判別確率:p=" << pr << "(誤判別確率:1-p=" << 1.0L - pr << ")\n";
  pr = ftop(f, var.df1(), var.df2());
  fout << "・HotellingのT2乗検定:T2乗=D2乗×" << var.n(0) << "×" << var.n(1) << "/" << var.nt() << "=" << t2 << "\n";
  fout << "  F=" << f << "  第1自由度=" << var.df1() << "  第2自由度=" << var.df2() << "  有意確率:p=" << pr << pind(pr) << endl;
}

void disanclass::put_dnsid(ofstream& fout, option& opt)
// DANSデータのタイトルとデータ識別情報を出力する
{
  idd tid(id.title(), 3);
  id.put_title(fout, formtag);
  id.put_stra(fout, formtag);
  if (opt.dns() == dnscoef)     // 判別係数出力
  {
    tid.assign_item(0, "項目名");
    if (opt.fin() > 0.0 || opt.pin() < 1.0) tid.assign_item(1, "判別係数(変数選択法)");
    else tid.assign_item(1, "判別係数(全変数選択)");
    tid.assign_item(2, "欠損補完用平均値");
    tid.put_item(fout, formtag);
  }
  else                          // 判別スコア出力
  {
    tid.assign_item(0, id.item(0));
    if (opt.fin() > 0.0 || opt.pin() < 1.0)
    {
      tid.assign_item(1, "判別スコア(変数選択法)");
      tid.assign_item(2, "群1に属す確率(変数選択法)");
    }
    else
    {
      tid.assign_item(1, "判別スコア(全変数選択)");
      tid.assign_item(2, "群1に属す確率(全変数選択)");
    }
    tid.put_item(fout, formtag);
  }
}

void disanclass::put_coef(ofstream& fout, const varclass& var)
// 判別スコア計算用判別係数をデータファイル形式で出力する
{
  datac cdata(3);
  cdata.assign(0, "切片");
  cdata.assign(1, var.a(0));
  cdata.assign(2, 0.0);
  cdata.put(fout);
  for (int i = 1; i < var.np(); ++i)
  {
    if (var.check(i))
    {
      cdata.assign(0, id.item(i));
      cdata.assign(1, var.a(i));
      cdata.assign(2, var.mid(i));
      cdata.put(fout);
    }
  }
}

void disanclass::put_score(ofstream& fout, option& opt, const varclass& var)
// 判別スコアと群1に属す確率をデータファイル形式で出力する
{
  double pi1 = (double)(var.n(0)) / (double)(var.nt());
  double cf = (double)(var.n(1)) / (double)(var.n(0));
  if (opt.pi() != 0.0)
  { pi1 = opt.pi();  cf = (1.0 - opt.pi()) / opt.pi(); }
  datac cdata(3);
  for (auto& cdatai : datalst)
  {
    cdata.assign(0, cdatai.sdata(0));
    double z = BLANKf;
    int gidx = var.calscore(z, opt, cdatai);
    if (gidx < 0)
    { cdata.assign(1, BLANKs);  cdata.assign(2, BLANKs); }
    else
    { cdata.assign(1, z);  cdata.assign(2, 1.0 / (1.0 + cf * exp(-z))); }
    cdata.put(fout);
  }
}

void disanclass::plot_disngp(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, opt, var);
    plot_log(fout, opt, var);
    delml = true;
  }
  if (opt.ngp("lci"))           // ロジスティック曲線の信頼区間
  {
    if (delml)
    { fout << NGPdelm << "\n"; }
    else
    {
      fout << "-s1 -dline -w80 -c255,0,0 -zy-0.1,1.1,0.2\n";
      delml = true;
    }
    plot_ngpid(fout, opt, var);
    plot_lci(fout, opt, var);
  }
  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, opt, var);
    plot_score(fout, opt, var);
  }
}

void disanclass::plot_ngpid(ofstream& fout, option& opt, const varclass& var)
// Ngraph用データのタイトルとデータ識別情報を出力する
{
  fout << "%";  puttitle(fout, Title);
  fout << "%\n";
  id.put_title(fout, formngp);
  id.put_stra(fout, formngp);
  idd tid(id);
  tid.set(id.title(), id.ni() + (opt.ixlst()).size()); // 項目数を変更する
  set_itemix(tid, opt, id.ni());
  fout << "%判別項目:" << tid.item(0) << "\n";
  for (size_t i = 1; i < tid.ni(); ++i)
  {
    if (var.check(i))
    {
      fout << "%変数x" << setw(3) << i;
      fout << ":" << tid.item(i) << "\n";
    }
  }
  if (opt.pin() < 1.0) fout << "%取り込み基準:p値＜Pin=" << opt.pin() << "  追い出し基準:p値≧Pout=" << opt.pout() << "\n";
  else fout << "%取り込み基準:偏F値＞Fin=" << opt.fin() << "  追い出し基準:偏F値≦Fout=" << opt.fout() << "\n";
  if (opt.xno() > 0)          // X軸にする変数指定あり
  {
    fout << "%X軸にする変数:" << tid.item(opt.xno()) << "\n";
    for (size_t i = 1; i < tid.ni(); ++i)
    {
      if (var.check(i))
      {
        if (opt.x(i) != BLANKf) fout << "%X軸にしない変数:" << tid.item(i) << " 指定値=" << opt.x(i) << "\n";
      }
    }
  }
}

void disanclass::plot_log(ofstream& fout, option& opt, const varclass& var)
// Ngraph用ロジスティック曲線を出力する
{
  datac cdata(id.ni());
  cdata.assign(0, var.a(0));    // 切片
  for (size_t i = 1; i < id.ni(); ++i)
  {
    if (opt.x(i) == BLANKf) cdata.assign(i, var.mid(i));
    else cdata.assign(i, opt.x(i));
  }
  double pi1 = (double)(var.n(0)) / (double)(var.n(0) + var.n(1));
  double cf = (double)(var.n(1)) / (double)(var.n(0));
  if (opt.pi() != 0.0) { pi1 = opt.pi();  cf = (1.0 - opt.pi()) / opt.pi(); }
  fout << "%ロジスティック曲線\n";
  fout << "%群1の事前確率=" << pi1 << "\n";
  double xmin = 0.0, xmax = 0.0;
  double xrange = var.calrange(xmin, xmax, opt, cdata);
  double dx = xrange / (double)RDdef;
  cdata.assign(0, opt.g(0));
  double x = xmin;
  for (int i = 0; i <= RDdef; ++i, x += dx)
  {
    double z = BLANKf;
    if (opt.xno() < 1)
    { z = x; }
    else
    {
      cdata.assign(opt.xno(), x);
      var.calscore(z, opt, cdata);
    }
    double pr = 1.0 / (1.0 + cf * exp(-z));
    fout << x << " " << pr << "\n";
  }
}

void disanclass::plot_lci(ofstream& fout, option& opt, const varclass& var)
// Ngraph用ロジスティック曲線の信頼区間を出力する
{
  datac cdata(id.ni());
  cdata.assign(0, var.a(0));    // 切片
  for (size_t i = 1; i < id.ni(); ++i)
  {
    if (opt.x(i) == BLANKf) cdata.assign(i, var.mid(i));
    else cdata.assign(i, opt.x(i));
  }
  double pi1 = (double)(var.n(0)) / (double)(var.nt());
  double cf = (double)(var.n(1)) / (double)(var.n(0));
  if (opt.pi() != 0.0)
  { pi1 = opt.pi();  cf = (1.0 - opt.pi()) / opt.pi(); }
  fout << "%ロジスティック曲線の" << opt.cc() * 100.0 << "%信頼区間\n";
  fout << "%群1の事前確率=" << pi1 << "\n";
  cdata.assign(0, opt.g(0));
  double t = ptoz(1.0 - opt.cc());
  double ci = t * sqrt((1.0 / (double)(var.n(0)) + 1.0 / (double)(var.n(1))) * var.d2());
  fout << "%信頼区間下限\n";
  double xmin = 0.0, xmax = 0.0;
  double xrange = var.calrange(xmin, xmax, opt, cdata);
  double dx = xrange / (double)RDdef;
  double x = xmin;
  for (int i = 0; i <= RDdef; ++i, x += dx) // 信頼区間下限データ出力
  {
    double z = BLANKf;
    if (opt.xno() < 1)
    { z = x; }
    else
    {
      cdata.assign(opt.xno(), x);
      var.calscore(z, opt, cdata);
    }
    double pr = 1.0 / (1.0 + cf * exp(-z + ci));
    fout << x << " " << pr << "\n";
  }
  fout << NGPdelm2 << "\n";
  fout << "%信頼区間上限\n";
  x = xmin;
  for (int i = 0; i <= RDdef; ++i, x += dx) // 信頼区間上限データ出力
  {
    double z = BLANKf;
    if (opt.xno() < 1)
    { z = x; }
    else
    {
      cdata.assign(opt.xno(), x);
      var. calscore(z, opt, cdata);
    }
    double pr = 1.0 / (1.0 + cf * exp(-z - ci));
    fout << x << " " << pr << "\n";
  }
}

void disanclass::plot_score(ofstream& fout, option& opt, const varclass& var)
// Ngraph用判別スコアを出力する
{
  fout << "%判別スコア：" << id.item(0) << "=";
  if (opt.g(0).empty()) fout << opt.g(1) << " 以外";
  else fout << opt.g(0);
  fout << "\n";
  for (auto& cdatai : datalst)
  {
    double x = BLANKf, z = BLANKf;
    int gidx = var.calscore(z, opt, cdatai);
    if (gidx == 0)
    {
      if (opt.xno() < 1) x = z;
      else x = cdatai.fdata(opt.xno());
      if (x != BLANKf) fout << x << " 1\n";
    }
  }
  fout << NGPdelm2 << "\n";
  fout << "%判別スコア：" << id.item(0) << "=";
  if (opt.g(1).empty()) fout << opt.g(0) << " 以外";
  else fout << opt.g(1);
  fout << "\n";
  for (auto& cdatai : datalst)
  {
    double x = BLANKf, z = BLANKf;
    int gidx = var.calscore(z, opt, cdatai);
    if (gidx == 1)
    {
      if (opt.xno() < 1) x = z;
      else x = cdatai.fdata(opt.xno());
      if (x != BLANKf) fout << x << " 0\n";
    }
  }
}

int disanclass::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 disanclass::put(option& opt) // 集計結果を出力する
{
  varclass var(id.ni(), opt);   // 変数オブジェクト
  var.cal(opt, datalst);        // 各種の計算実行
  if (opt.res("res"))
  {
    put_resid(cout, opt);
    if (opt.res("statg")) put_statg(cout, opt, var);
    if (opt.res("stata"))
    {
      cout << "\n○全体(SD・SE:群内)";
      putmsdvec(cout, var.nt(), var.mevec(), var.vrvec());
      putmat(cout, "・群内相関行列", matr, var.nt() - 1, var.dmat(), var.nxix());
    }
    if (var.resvar() == reserr)
    { cout << "\n※判別関数計算不可能！\n\n" << endl; }
  }
  if (var.resvar() == reserr) return;
  var.stepsel(opt.fin(), opt.fout(), opt.pin(), opt.pout());
  if (opt.res("res") && opt.res("selx")) put_distab(cout, opt, opt.fin(), opt.fout(), opt.pin(), opt.pout(), var);
  if (opt.dnsname() != "")
  {
    ofstream fout;
    if (ofopen(fout, (opt.dnsname()).c_str()) == false) return; // ファイルオープン失敗
    put_dnsid(fout, opt);
    if (opt.dns() == dnscoef) put_coef(fout, var);
    else put_score(fout, opt, var);
  }
  if (opt.ngpname() != "") plot_disngp(opt, var);
  if (opt.res("res"))
  {
    if (opt.res("allx") && (opt.res("selx") == false || var.df1() < var.nxix()))
    {
      if (var.df1() < var.np()) var.stepsel(0.0, 0.0, 1.0, 1.0);
      put_distab(cout, opt, 0.0, 0.0, 1.0, 1.0, var);
    }
    cout << "\n" << endl;
  }
}

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