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

constexpr double CCmra = 0.5;   // 省略時信頼係数(CCdefとは別の値)

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;               // プログラム名
  double v_cc;                  // 信頼係数
  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; }
  double cc() const             // 信頼係数を返す
  { return v_cc; }
  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(), v_cc(CCmra), 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); // プログラム名削除
}

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() ? CCmra : strtodf(opt[1]);
    if (v_cc < 0.0 || 1.0 <= v_cc) err = badopt(pname(), optstr, USAGE, 0.0, 1.0);
  }
  else if (opt[0] == "f")       // オプション指定ファイル
  { err = fgetopt(*this, opt[1], true); }
  else if (opt[0] == "h")       // オンラインヘルプ表示
  { err = help(pname(), USAGE, opt[1]); }
  else if (opt[0] == "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_ni;                     // 項目数
  ldvec v_cvdveclst;            // 交差検証法の検証結果リスト
  restype v_resvar;             // 処理結果の種類
  restype cal_mra(int& nt, mat<double>& smat, double& syy, lditr& lditrlstt) const;
  // 重回帰分析の計算をする
  int cv_mra(option& opt, const mat<double>& smat, int dfr, lditr& lditrlstv) const;
  // 重回帰分析結果でデータリストを検証してyが予測限界内に入った例数を返す
  void set_cvdlst(int vno, int ny, int yfit, const mat<double>& smat, double syy);
  // 検証結果を検証結果リストに追加して要素数を返す
  void cal_cv(option& opt, lditr& lditrlst, int r, int k);
  // 重回帰分析の交差検証を行う
  double cal_yest(datac& cdata, const mat<double>& smat) const;
  // y推定値を計算して返す
  double cal_pi(datac& cdata, const mat<double>& smat, double vr, double tcc) const;
  // yの予測限界幅を計算して返す
public:
  varclass(size_t ni) : v_ni(ni), v_cvdveclst(), v_resvar(resnon)
  // コンストラクタ
  {}
  ~varclass()                   // デストラクタ
  {}
  int ni() const                // 項目数を返す
  { return v_ni; }
  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 = ni() + 1, nt = 0;
  mat<double> smat(dim);        // 積和行列
  double syy = 0.0;
  v_resvar = cal_mra(nt, smat, syy, lditrlst);
  if (v_resvar == reserr) return;
  int yfit0 = cv_mra(opt, smat, nt - ni(), lditrlst); //yが予測限界内に入った例数 
  set_cvdlst(0, nt, yfit0, smat, syy);
  int r = opt.r();              // 検証群のデータ数
  if (r < 1 || nt - r < ni())    // 交差検証実行せず
  {
    v_resvar = reserr;
    return;
  }
  int k = nt / r;               // 検証群数
  if (nt % r > 0) ++k;          // 切り上げる
  cal_cv(opt, lditrlst, r, k);
  v_resvar = resok;
}

restype varclass::cal_mra(int& nt, mat<double>& smat, double& syy, lditr& lditrlstt) const
// 重回帰分析の計算をする
{
  int dim = ni() + 1;
  vecdbl fdatavec(dim, 0.0);
  fdatavec[0] = 1.0;
  for (lditr_itr lditr = lditrlstt.begin(); lditr != lditrlstt.end(); ++lditr)
  {
    for (int i = 1; i < ni(); ++i) fdatavec[i] = (*lditr)->fdata(i); // x
    fdatavec[ni()] = (*lditr)->fdata(0); // y
    for (int i = 0; i < dim; ++i)
    {
      double xyi = fdatavec[i];
      for (int j = i; j < dim; ++j) smat.add(i, j, xyi * fdatavec[j]);
    }
  }
  nt = (int)(smat(0, 0));
  if (nt < 1) return reserr;    // データなし
  for (int i = 0; i < ni(); ++i)
  {
    for (int j = i + 1; j < dim; ++j) smat.assign(j, i, smat(i, j));
  }
  smat.sweep(0);                // 単純積和行列→積和行列
  syy = smat(ni(), ni());
  for (int i = 1; i < ni(); ++i) smat.sweep(i); // 行列を掃き出す
  if (sgn(syy) <= 0) return reserr;
  return resok;
}

int varclass::cv_mra(option& opt, const mat<double>& smat, int dfr, lditr& lditrlstv) const
// 重回帰分析結果でデータリストを検証してyが予測限界内に入った例数を返す
{
  double vr = (dfr > 0) ? smat(ni(), ni()) / (double)dfr : 0.0;
  double tcc = ptot(1.0 - opt.cc(), dfr);
  int yfit = 0;
  for (lditr_itr lditr = lditrlstv.begin(); lditr != lditrlstv.end(); ++lditr)
  {
    double y = (*lditr)->fdata(0);
    double yest = cal_yest(*(*lditr), smat);
    double pi = cal_pi(*(*lditr), smat, vr, tcc);
    double ylow = yest - pi, yup = yest + pi;
    if (ylow <= y && y <= yup) ++yfit;
  }
  return yfit;                  // yが予測限界内に入った例数を返す
}

void varclass::set_cvdlst(int vno, int ny, int yfit, const mat<double>& smat, double syy)
// 検証結果を検証結果リストに追加して要素数を返す
{
  vecdbl cvdvec(ni() + 4, 0.0);
  cvdvec[0] = (double)((vno < 1) ? ni() + 1 : vno);
  for (int i = 0; i < ni(); ++i) cvdvec[i + 1] = smat(i, ni());
  cvdvec[ni() + 1] = 1.0 - smat(ni(), ni()) / syy;
  cvdvec[ni() + 2] = (double)ny;
  cvdvec[ni() + 3] = (double)yfit;
  v_cvdveclst.push_back(cvdvec);
}

void varclass::cal_cv(option& opt, lditr& lditrlst, int r, int k)
// 重回帰分析の交差検証を行う
{
  int dim = ni() + 1, nt = 0;
  mat<double> smat(dim, dim); // 積和行列
  double syy = 0.0;           // yの平方和
  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);
    }
    smat = 0.0;  syy = 0.0;
    if (cal_mra(nt, smat, syy, lditrlstt) != reserr)
    {
      int yfit = cv_mra(opt, smat, nt - ni(), lditrlstv);
      int ny = lditrlstv.size();
      set_cvdlst(i + 1, ny, yfit, smat, syy);
    }
    lditrlstt.clear();  lditrlstv.clear();
  }
}

double varclass::cal_yest(datac& cdata, const mat<double>& smat) const
// y推定値を計算して返す
{
  double yest = smat(0, ni());
  for (int i = 1; i < ni(); ++i) yest += smat(i, ni()) * cdata.fdata(i);
  return yest;                  // y推定値を返す
}

double varclass::cal_pi(datac& cdata, const mat<double>& smat, double vr, double tcc) const
// yの予測限界幅を計算して返す
{
  int dim = ni() + 1;
  vec<double> xv(dim);
  for (int i = 1; i < ni(); ++i) xv.assign(i, cdata.fdata(i));
  xv.assign(0, 1.0);  xv.assign(ni(), 0.0);
  vec<double> tv(xv.trn() * smat);
  double vy = vr * (1.0 + inp(tv, xv));
  return tcc * sqrt(vy);        // yの予測限界幅を返す
}

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

void mracvclass::put_cv(ostream& fout, option& opt, varclass& var)
// 重回帰分析の交差検証結果をデータファイル形式で出力する
{
  id.put_title(fout, formtag);
  id.put_stra(fout, formtag);
  size_t ni = id.ni();
  idd tid(id.title(), ni + 4);
  tid.assign_item(0, "パラメータ数/検証群番号");
  tid.assign_item(1, "b0:切片");
  for (size_t i = 1; i < ni; ++i)
  {
    ostringstream ostr;
    ostr << "b" << i << ":" << id.item(i);
    tid.assign_item(1 + i, ostr.str());
  }
  ostringstream ostr;
  ostr << id.item(0) << "の重寄与率:R2乗";
  tid.assign_item(ni + 1, ostr.str());
  tid.assign_item(ni + 2, "検証群の例数");
  ostr.str("");
  ostr << opt.cc() * 100.0 << "%予測限界内の例数";
  tid.assign_item(ni + 3, ostr.str());
  tid.put_item(fout, formtag);
  int prec0 = fout.precision(); // 現在の浮動小数点精度(有効数字)を保存する
  for (auto& cvi : var.cvdveclst())
  {
    fout << cvi[0];
    fout.precision(16);
    for (size_t i = 1; i <= ni + 1; ++i) fout << CHRdelm << cvi[i];
    fout.precision(prec0);
    for (size_t i = ni + 2; i < tid.ni(); ++i) fout << CHRdelm << cvi[i];
    fout << endl;
  }
}

int mracvclass::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("項目数が足りません！");
  nc = id.nc();                 // ケース数変更用
  ldata_itr ditr = datalst.begin();
  while (ditr != datalst.end())
  {
    bool blank = false;
    for (size_t j = 0; blank == false && j < id.ni(); ++j)
    {
      if (ditr->fdata(j) == BLANKf) blank = true;
    }
    if (blank)                  // ケースを削除する
    { ditr = datalst.erase(ditr);  --nc; }
    else
    { ++ditr; }
  }
  if (nc < (int)(id.nc())) id.nc(nc);  // ケース数を減らす
  return nc;                           // 有効なケース数を返す
}

void mracvclass::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, opt, var);
}

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