// ldacv.C Version 7.43  Copyleft (c) 杉本典夫
// Time-stamp: <2025-10-31 09:10:07 tac>
//
// 最初の項目を判別項目、以後の項目を説明変数とした判別分析の交差検証を行う。
//
// 使用ライブラリ:libdans.so, libfun.so
// コンパイル法  :g++ -O -s -W ldacv ldacv.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 =            // プログラム使用法表示用
{
  "ldacv  Version 7.43          Copyleft(Vo\\oV) 2025 杉本典夫",
  "書式:ldacv [<オプション> …] [<データファイル> …]",
  "機能:最初の項目を判別項目、以後の項目を説明変数とした",
  "     判別分析の交差検証を行います。",
  "<オプション>",
  "  -f=<ファイル>:指定したファイルからオプションを読み込みます。",
  "  -g{1|2}[=<群コード>]:群1または群2のコードを指定します。",
  "  -h[=<コマンド名>]:オンラインマニュアルを表示します。(-h:-h=man)",
  "  -r=<データ数>:検証群のデータ数を指定します。"
};

using lditr = list<ldata_itr>; // データリスト反復子リストの別名定義
using lditr_itr = list<ldata_itr>::iterator; // 反復子リスト反復子の別名定義
using ldvec = list<vecdbl>;     // ベクトルデータの別名定義
using ldvec_itr = list<vecdbl>::iterator; // 反復子の別名定義

class option                    // オプションクラス
{
private:
  string v_path;                // プログラムのパス
  string v_pname;               // プログラム名
  vecstr gvec;                  // 群コード
  int v_r;                      // 検証群のデータ数
public:
  option(const char* c_pname = nullptr); // コンストラクタ
  ~option()                              // デストラクタ
  {}
  const string& path() const    // プログラムのパスを返す
  { return v_path; }
  const string& pname() const   // プログラム名を返す
  { return v_pname; }
  const string& g(int idx) const // 群コードを返す
  { return gvec[idx]; }
  int r() const                 // 検証群データ数を返す
  { return v_r; }
  void setdef()                 // 省略時オプション指定を設定する
  { getoptdef(*this); }
  void set(const string& optstr, char devc = OPTdev);
  // オプション指定によってオプションを設定する
};

option::option(const char* c_pname) : v_path(), v_pname(), gvec(2), v_r(1)
// コンストラクタ
{
  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); // プログラム名削除
  gvec[0] = "1";
}

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] == "f")            // オプション指定ファイル
  { err = fgetopt(*this, opt[1], true); }
  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] == "r")       // 検証群データ数の指定
  {
    char* ptr = nullptr;
    if (opt[1].empty()) v_r = 1;
    else v_r = int(strtol(opt[1].c_str(), &ptr, 10));
    if (v_r < 0) err = badopt(pname(), optstr, USAGE, 1.0, (double)INT_MAX);
  }
  else                          // 間違ったオプション指定
  { err = nonopt(pname(), optstr, USAGE); }
  if (err > 0) throw help_except("ヘルプ表示");
  if (err < 0) throw option_error(string("オプション指定違反→") + optstr);
}

class varclass                  // 変数クラス
{
private:
  int v_np;                     // パラメータ数
  ldvec v_cvdveclst;            // 交差検証法の結果リスト
  restype v_resvar;             // 計算結果の種類
  restype cal_lda(int& nt, mat<double>& dmat, lditr& lditrlstt);
  // 判別分析の計算をする
  int cv_lda(int& nd, int& tp, int& nn, int& fn, const mat<double>& dmat, lditr& lditrlstv) const;
  // 判別分析結果でデータリストを検証して検証例数を返す
  int dis_score(double& z, datac& cdata, const mat<double>& dmat) const;
  // 判別スコアを計算し、群インデックスを返す
  void set_cvdlst(int vno, int nd, int tp, int nn, int fn, const mat<double>& dmat);
  // 検証結果を検証結果リストに追加する
  void cal_cv(mat<double>& dmat, lditr& lditrlst, int r, int k);
  // 判別分析の交差検証を行う
public:
  varclass(size_t ni) : v_np(ni), v_cvdveclst(), v_resvar(resnon)
  // コンストラクタ
  {}
  ~varclass()                   // デストラクタ
  {}
  int np() const                // パラメータ数を返す
  { return v_np; }
  ldvec& cvdveclst()            // 交差検証法の結果リストを返す
  { return v_cvdveclst; }
  restype resvar() const        // 計算結果の種類を返す
  { return v_resvar; }
  void cal(option& opt, lditr& lditrlst); // データを集計して各種統計量を求める
};

void varclass::cal(option& opt, lditr& lditrlst)
// データを集計して各種統計量を求める
{
  int dim = np() + 1, nt = 0;
  mat<double> dmat(dim);        // 群合計積和行列
  v_resvar = cal_lda(nt, dmat, lditrlst);
  if (v_resvar == reserr) return;
  int r = opt.r();              // 検証群のデータ数
  if (r < 1 || nt - r < np())   // 交差検証実行せず
  {
    v_resvar = reserr;
    return;
  }
  int k = nt / r;               // 検証群数
  if (nt % r > 0) ++k;          // 切り上げる
  int nd0 = 0, tp0 = 0, nn0 = 0, fn0 = 0;
  cv_lda(nd0, tp0, nn0, fn0, dmat, lditrlst);
  set_cvdlst(0, nd0, tp0, nn0, fn0, dmat);
  cal_cv(dmat, lditrlst, r, k);
  v_resvar = resok;
}

restype varclass::cal_lda(int& nt, mat<double>& dmat, lditr& lditrlstt)
// 判別分析の計算をする
{
  int dim = np() + 1;
  vector<mat<double>> smatvec(2); // 積和行列ベクトル
  smatvec[0].resize(dim, dim);  smatvec[1].resize(dim, dim);
  smatvec[0] = 0.0;  smatvec[1] = 0.0;
  vecint nvec(2, 0);            // 例数ベクトル
  vecdbl fdatavec(np(), 0.0);
  fdatavec[0] = 1.0;
  for (lditr_itr lditr = lditrlstt.begin(); lditr != lditrlstt.end(); ++lditr)
  {
    int idx = (int)((*lditr)->fdata(0));
    for (int i = 1; i < np(); ++i) fdatavec[i] = (*lditr)->fdata(i);
    for (int i = 0; i < np(); ++i)
    {
      double x = fdatavec[i];
      for (int j = i; j < np(); ++j) smatvec[idx].add(i, j, x * fdatavec[j]);
    }
  }
  nvec[0] = (int)(smatvec[0](0, 0));  nvec[1] = (int)(smatvec[1](0, 0));
  nt = nvec[0] + nvec[1];
  if (nt < 1) return reserr;   // データなし
  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);
    }
  }
  if (nvec[0] < 1 || nvec[1] < 1 || nt < np()) return reserr;
  dmat = smatvec[0] + smatvec[1];
  int df = nt - 2;
  dmat *= 1.0 / (double)df;
  dmat.assign(0, 0, 0.0);  dmat.assign(0, np(), 0.0);
  dmat.assign(np(), 0, 0.0);  dmat.assign(np(), np(), 0.0);
  for (int i = 1; i < np(); ++i)
  {
    dmat.assign(0, i, (smatvec[0](0, i) + smatvec[1](0, i)) * 0.5);
    dmat.assign(i, 0, -dmat(i, 0));
    dmat.assign(i, np(), smatvec[0](0, i) - smatvec[1](0, i));
    dmat.assign(np(), i, dmat(i, np()));
  }
  for (int i = 1; i < np(); ++i) dmat.sweep(i);
  return resok;
}

int varclass::cv_lda(int& nd, int& tp, int& nn, int& fn, const mat<double>& dmat, lditr& lditrlstv) const
// 判別分析結果でデータリストを検証して検証例数を返す
{
  nd = tp = nn = fn = 0;
  for (lditr_itr lditr = lditrlstv.begin(); lditr != lditrlstv.end(); ++lditr)
  {
    double z = 0.0;
    int gidx = dis_score(z, *(*lditr), dmat);
    if (gidx == 0)              // 疾患群
    {
      ++nd;
      if (z >= 0.0) ++tp;
    }
    else if (gidx == 1)         // 正常群
    {
      ++nn;
      if (z < 0.0) ++fn;
    }
  }
  return nd + nn;               // 合計例数を返す
}

int varclass::dis_score(double& z, datac& cdata, const mat<double>& dmat) const
// 判別スコアを計算し、群インデックスを返す
{
  z = dmat(0, np());
  for (int i = 1; i < np(); ++i) z += dmat(i, np()) * cdata.fdata(i);
  return (int)(cdata.fdata(0)); // 群インデックスを返す
}

void varclass::set_cvdlst(int vno, int nd, int tp, int nn, int fn, const mat<double>& dmat)
// 検証結果を検証結果リストに追加する
{
  vecdbl cvdvec(np() + 8, 0.0);
  cvdvec[0] = (double)((vno < 1) ? np() + 1 : vno);
  for (int i = 0; i < np(); ++i) cvdvec[1 + i] = dmat(i, np());
  cvdvec[np() + 1] = -dmat(np(), np());
  cvdvec[np() + 2] = (double)nd;
  cvdvec[np() + 3] = (double)tp;
  cvdvec[np() + 4] = (double)nn;
  cvdvec[np() + 5] = (double)fn;
  cvdvec[np() + 6] = (double)(nd + nn);
  cvdvec[np() + 7] = (double)(tp + fn);
  v_cvdveclst.push_back(cvdvec);
}

void varclass::cal_cv(mat<double>& dmat, lditr& lditrlst, int r, int k)
// 判別分析の交差検証を行う
{
  lditr lditrlstt, lditrlstv; // トレーニングデータとバリデーションデータ反復子
  for (int i = 0; i < k; ++i)
  {
    int sidx = r * i, eidx = r * (i + 1) - 1;
    int j = 0;
    for (lditr_itr lditr = lditrlst.begin(); lditr != lditrlst.end(); ++lditr, ++j)
    {
      if (sidx <= j && j <= eidx) lditrlstv.push_back(*lditr);
      else lditrlstt.push_back(*lditr);
    }
    dmat = 0.0;
    int nt = 0;
    if (cal_lda(nt, dmat, lditrlstt) > 0)
    {
      int nd = 0, tp = 0, nn = 0, fn = 0;
      cv_lda(nd, tp, nn, fn, dmat, lditrlstv);
      set_cvdlst(i + 1, nd, tp, nn, fn, dmat);
    }
    lditrlstt.clear();  lditrlstv.clear();
  }
}

class ldacvclass                // ldacvメインクラス
{
private:
  string fname;                 // ファイル名
  idd id;                       // データ識別情報
  ldata datalst;                // データリスト
  int search_g(option& opt, const string& gcode) const;
  // 群コードを検索する
  void put_cv(ostream& fout, varclass& var);
  // 判別分析の交差検証結果をデータファイル形式で出力する
public:
  ldacvclass() : fname(), id(), datalst() // コンストラクタ
  {}
  ~ldacvclass()                 // デストラクタ
  {}
  int get(const char* c_fname, option& opt); // データを読み込み累積する
  void put(option& opt);                     // 集計結果を出力する
};

int ldacvclass::search_g(option& opt, const string& gcode) const
// 群コードを検索する
{
  if (gcode == BLANKs) return -1; // 欠損データはヒットさせない
  int idx = -1;
  if ((opt.g(0)).empty())       // 群1のコード指定がない場合
  { idx = (gcode == opt.g(1)) ? 1 : 0;  }
  else if (gcode == opt.g(0))   // 群1のコード
  { idx = 0; }
  else if ((opt.g(1)).empty())  // 群2のコード指定がない場合
  { idx = 1; }
  else if (gcode == opt.g(1))   // 群2のコード
  { idx = 1; }
  return idx;                   // -1:ヒットせず,0:群1,1:群2
}

void ldacvclass::put_cv(ostream& fout, varclass& var)
// 判別分析の交差検証結果をデータファイル形式で出力する
{
  id.put_title(fout, formtag);
  id.put_stra(fout, formtag);
  size_t ni = id.ni();
  idd tid(id.title(), ni + 8);
  tid.assign_item(0, "パラメータ数/検証群番号");
  tid.assign_item(1, "a0:切片");
  for (size_t i = 1; i < ni; ++i)
  {
    ostringstream ostr;
    ostr << "a" << i << ":" << id.item(i);
    tid.assign_item(1 + i, ostr.str());
  }
  tid.assign_item(ni + 1, "Mahalanobisの汎距離:D2乗");
  tid.assign_item(ni + 2, "群1(疾患群)の例数");
  tid.assign_item(ni + 3, "群1(疾患群)の正判別例数(TP)");
  tid.assign_item(ni + 4, "群2(正常群)の例数");
  tid.assign_item(ni + 5, "群2(正常群)の正判別例数(FN)");
  tid.assign_item(ni + 6, "2群合計例数");
  tid.assign_item(ni + 7, "2群合計正判別例数(AC)");
  tid.put_item(fout, formtag);
  int prec0 = fout.precision(); // 現在の浮動小数点精度(有効数字)を保存する
  for (ldvec_itr cvitr = (var.cvdveclst()).begin(); cvitr != (var.cvdveclst()).end(); ++cvitr)
  {
    fout << (*cvitr)[0];
    fout.precision(16);
    for (size_t i = 1; i <= ni + 1; ++i) fout << CHRdelm << (*cvitr)[i];
    fout.precision(prec0);
    for (size_t i = ni + 2; i < tid.ni(); ++i) fout << CHRdelm << (*cvitr)[i];
    fout << endl;
  }
}

int ldacvclass::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.g(0).empty() && opt.g(1).empty()) throw runtime_error("群コードが指定されていません！");
  if (id.ni() < 2) throw runtime_error("項目数が足りません！");
  nc = (int)(id.nc());          // ケース数変更用
  ldata_itr ditr = datalst.begin();
  while (ditr != datalst.end())
  {
    bool blank = false;
    for (size_t j = 1; blank == false && j < id.ni(); ++j)
    {
      if (ditr->fdata(j) == BLANKf) blank = true;
    }
    if (ditr->sdata(0) == BLANKs || blank) // 欠損のためケースを削除する
    { ditr = datalst.erase(ditr);  --nc; } // 削除した次の要素の位置を返す
    else
    {
      int idx = search_g(opt, ditr->sdata(0));
      if (idx >= 0)             // 群コードをインデックスに変える
      {
        ditr->assign(0, (double)idx);
        ++ditr;
      }
      else                      // ケースを削除する
      { ditr = datalst.erase(ditr);  --nc; }
    }
  }
  if (nc < (int)(id.nc())) id.nc(nc);  // ケース数を減らす
  return nc;                    // 有効なケース数を返す
}

void ldacvclass::put(option& opt) // 集計結果を出力する
{
  lditr lditrlst;               // データリスト反復子リスト
  for (ldata_itr ditr = datalst.begin(); ditr != datalst.end(); ++ditr) lditrlst.push_back(ditr);
  varclass var(id.ni());        // 変数オブジェクト
  var.cal(opt, lditrlst);
  if (var.resvar() == reserr) return;
  put_cv(cout, var);
}

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