// cvjack.C Version 7.42  Copyleft (c) 杉本典夫
// Time-stamp: <2025-11-03 10:58:18 tac>
//
// 交差検証データのジャックナイフ推定値を計算する。
//
// 使用ライブラリ:libdans.so, libfun.so
// コンパイル法  :g++ -O -s -W cvjack cvjack.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 =            // プログラム使用法表示用
{
  "cvjack  Version 7.42         Copyleft(Vo\\oV) 2025 杉本典夫",
  "書式:cvjack [<オプション> …] [<データファイル> …]",
  "機能:交差検証データのジャックナイフ推定値を計算します。",
  "<オプション>",
  "  -cc[=<信頼係数>]:信頼係数を指定します。(0≦信頼係数＜1)",
  "  -f=<ファイル>:指定したファイルからオプションを読み込みます。",
  "  -h[=<コマンド名>]:オンラインマニュアルを表示します。(-h:-h=man)",
  "  -r=<検証群のデータ数>:検証群のデータ数を指定します。"
};

const string Title = "=== 交差検証(cross-validation) ==="; // 表題

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);
  // オプション指定によってオプションを設定する
  int check(size_t ni);         // オプションのチェックとセット
};

option::option(const char* c_pname) : v_path(), v_pname(), v_cc(CCdef), 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() ? 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] == "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 < 1) 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);
}

int option::check(size_t ni)    // オプションのチェックとセット
{
  if (ni < 2) return disperr(pname(), "項目数が足りません！");
  return 0;                     // エラーなし
}

class jackest                   // ジャックナイフ推定値クラス
{
private:
  int v_k;                      // サブグループ数
  double v_h0;                  // 母数θ0
  double v_sigh;                // Σθi
  double v_sighm;               // Σθmi
  double v_sighm2;              // Σθmi2乗
public:
  explicit jackest() : v_k(0), v_h0(0.0), v_sigh(0.0), v_sighm(0.0), v_sighm2(0.0)
  // デフォールトコンストラクタ
  {}
  jackest(int k0, double h00 = 0.0) : v_k(k0), v_h0(h00), v_sigh(0.0), v_sighm(0.0), v_sighm2(0.0)
  // コンストラクタ
  {}
  jackest(const jackest& jackest0) : v_k(jackest0.v_k), v_h0(jackest0.v_h0), v_sigh(jackest0.v_sigh), v_sighm(jackest0.v_sighm), v_sighm2(jackest0.v_sighm2)
  // コピーコンストラクタ
  {}
  virtual ~jackest()            // デストラクタ
  {}
  const jackest& operator=(const jackest& jackest0); // コピー代入演算子
  int k() const                 // サブグループ数を返す
  { return v_k; }
  int k(int k0)                 // サブグループ数を設定する
  { return (v_k = k0); }
  double h0() const             // 母数θ0を返す
  { return v_h0; }
  double h0(double h00)         // 母数θ0を設定する
  { return (v_h0 = h00); }
  double sigh() const           // Σθiを返す
  { return v_sigh; }
  double sighm() const          // Σθmiを返す
  { return v_sighm; }
  double sighm2() const         // Σθmi2乗を返す
  { return v_sighm2; }
  void sigma(double hi);        // 母数の累積集計
  double h() const              // 母数の推定値を返す
  { return h0() * (double)(k()) - sigh() * (double)(k() - 1) / (double)(k()); }
  double vh() const;            // 推定値の分散を返す
  int df() const                // 自由度を返す
  { return (k() < 1) ? 0 : k() - 1; }
};

const jackest& jackest::operator=(const jackest& jackest0)
// コピー代入演算子
{
  if (&jackest0 == this) return *this; // 自分自身への代入を防ぐ
  v_k = jackest0.v_k;  v_h0 = jackest0.v_h0;
  v_sigh = jackest0.v_sigh;  v_sighm = jackest0.v_sighm;
  v_sighm2 = jackest0.v_sighm2;
  return *this;
}

void jackest::sigma(double hi)  // 母数の累積集計
{
  v_sigh += hi;
  double hmi = (double)(k()) * h0() - (double)(k() - 1) * hi;
  v_sighm += hmi;  v_sighm2 += hmi * hmi;
}

double jackest::vh() const      // 推定値の分散を返す
{  
  if (k() < 2) return 0.0;
  double v = (sighm2() - sighm() * sighm() / (double)k()) / (double)k() / (double)(k() - 1);
  const double Veps = FDeps * fmax(5.0, sighm2());
  return (v < Veps) ? 0.0 : v;
}

class varclass                  // 変数クラス
{
private:
  size_t v_np;                  // パラメータ数
  size_t v_nd;                  // 予測数
  vector<jackest> hjkvec;       // ジャックナイフ推定値ベクトル
  vecint n0vec;                 // 全体の合計例数ベクトル
  vecint nt0vec;                // 全体の正判別例数ベクトル
  vecint nvec;                  // 検証結果の合計例数ベクトル
  vecint ntvec;                 // 検証結果の正判別例数ベクトル
public:
  varclass(int np, int nd) : v_np(np), v_nd(nd), hjkvec(np), n0vec(nd, 0), nt0vec(nd, 0), nvec(nd, 0), ntvec(nd, 0)
  // コンストラクタ
  {}
  ~varclass()                   // デストラクタ
  {}
  size_t np() const             // パラメータを返す
  { return v_np; }
  size_t nd() const             // 予測数を返す
  { return v_nd; }
  void hjksigma(size_t idx, double hi) // 母数を累積する
  { if (idx < np()) hjkvec[idx].sigma(hi); }
  int hjkk(size_t idx) const    // サブグループ数を返す
  { return (idx < np()) ? hjkvec[idx].k() : 0; }
  int hjkk(size_t idx, int k0)  // サブグループ数をセットする
  { return (idx < np()) ? hjkvec[idx].k(k0) : 0.0; }
  double hjkh0(size_t idx) const // 母数θ0を返す
  { return (idx < np()) ? hjkvec[idx].h0() : 0.0; }
  double hjkh0(size_t idx, double h00) // 母数θ0をセットする
  { return (idx < np()) ? hjkvec[idx].h0(h00) : 0.0; }
  double hjkh(size_t idx) const // 母数推定値を返す
  { return (idx < np()) ? hjkvec[idx].h() : 0.0; }
  double hjkvh(size_t idx) const // 推定値の分散を返す
  { return (idx < np()) ? hjkvec[idx].vh() : 0.0; }
  int hjkdf(size_t idx) const   // 自由度を返す
  { return (idx < np()) ? hjkvec[idx].df() : 0; }
  int n0(size_t idx) const      // 全体の合計例数を返す
  { return (idx < nd()) ? n0vec[idx] : 0; }
  int n0(size_t idx, int n00)   // 全体の合計例数をセットする
  { return (idx < nd()) ? (n0vec[idx] = n00) : 0; }
  int nt0(size_t idx) const     // 全体の正判別例数を返す
  { return (idx < nd()) ? nt0vec[idx] : 0; }
  int nt0(size_t idx, int nt00) // 全体の正判別例数をセットする
  { return (idx < nd()) ? (nt0vec[idx] = nt00) : 0; }
  int n(size_t idx) const       // 検証結果の合計例数を返す
  { return (idx < nd()) ? nvec[idx] : 0; }
  int n(size_t idx, int n0)     // 検証結果の合計例数をセットする
  { return (idx < nd()) ? (nvec[idx] = n0) : 0; }
  int nt(size_t idx) const      // 検証結果の正判別例数を返す
  { return (idx < nd()) ? ntvec[idx] : 0; }
  int nt(size_t idx, int nt0)   // 検証結果の正判別例数をセットする
  { return (idx < nd()) ? (ntvec[idx] = nt0) : 0; }
};

class cvjackclass               // cvjackメインクラス
{
private:
  string fname;                 // ファイル名
  idd id;                       // データ識別情報
  ldata datalst;                // データリスト
  void cal_jack(varclass& var);
  // ジャックナイフ推定値を計算する
  void put_jack(ostream& fout, const option& opt, varclass& var);
  // ジャックナイフ推定値を出力する
public:
  cvjackclass() : fname(), id(), datalst() // コンストラクタ
  {}
  ~cvjackclass()                // デストラクタ
  {}
  int get(const char* c_fname, option& opt); // データを読み込み累積する
  void put(option& opt);                     // 集計結果を出力する
};

void cvjackclass::cal_jack(varclass& var)
// ジャックナイフ推定値を計算する
{
  ldata_itr dpos0 = datalst.begin(); // 最初のデータは全体の検証結果
  int k = id.nc() - 1;             // サブグループ数
  for (size_t i = 0; i < var.np(); ++i)
  {
    var.hjkk(i, k);
    var.hjkh0(i, dpos0->fdata(i + 1));
  }
  for (size_t i = 0; i < var.nd(); ++i)
  {
    size_t j = var.np() + i * 2;
    var.n0(i, (int)(dpos0->fdata(j + 1)));
    var.nt0(i, (int)(dpos0->fdata(j + 2)));
  }
  ++dpos0;
  for (ldata_itr dpos = dpos0; dpos != datalst.end(); ++dpos)
  {
    for (size_t i = 0; i < var.np(); ++i) var.hjksigma(i, dpos->fdata(i + 1));
    for (size_t i = 0; i < var.nd(); ++i)
    {
      size_t j = var.np() + i * 2;
      var.n(i, var.n(i) + (int)(dpos->fdata(j + 1)));
      var.nt(i, var.nt(i) + (int)(dpos->fdata(j + 2)));
    }
  }
}

void cvjackclass::put_jack(ostream& fout, const option& opt, varclass& var)
// ジャックナイフ推定値を出力する
{
  puttitle(fout, Title);
  fout << "\n";
  id.put_title(fout);
  id.put_stra(fout);
  fout << "・全体の検証結果\n";
  for (size_t i = 0; i < var.nd(); ++i)
  {
    size_t j = var.np() + i * 2;
    double pro = (var.n0(i) > 0) ? (double)(var.nt0(i)) / (double)(var.n0(i)) * 100.0 : 0.0;
    fout << "  " << id.item(j + 1) << "=" << var.n0(i);
    fout << "  " << id.item(j + 2) << "=" << var.nt0(i);
    fout << "(的中率=" << pro << "%)\n";
  }
  if (opt.r() > 1) fout << "\n・K-分割交差検証:検証群のデータ数=" << opt.r() << "　";
  else fout << "\n・1個抜き交差検証(LOOCV):";
  fout << "検証回数=" << var.hjkk(0) << "\n";
  for (size_t i = 0; i < var.nd(); ++i)
  {
    size_t j = var.np() + i * 2;
    double pro = (var.n(i) > 0) ? (double)(var.nt(i)) / (double)(var.n(i)) * 100.0 : 0.0;
    fout << "  " << id.item(j + 1) << "=" << var.n(i);
    fout << "  " << id.item(j + 2) << "=" << var.nt(i);
    fout << "(的中率=" << pro << "%)\n";
  }
  fout << "\n・ジャックナイフ推定値(jackknife estimate)\n";
  double t = ptot(1.0 - opt.cc(), var.hjkdf(0));
  for (size_t i = 0; i < var.np(); ++i)
  {
    double me = var.hjkh(i);
    double sei = (sgn(var.hjkvh(i)) > 0) ? sqrt(var.hjkvh(i)) : 0.0;
    double ci = t * sei;
    fout << "  " << id.item(i + 1) << "=" << me;
    double ml = me - ci, mu = me + ci;
    fout << " (" << opt.cc() * 100.0 << "%CI:" << ml << "-" << mu << ")" << "\n";
  }
}

int cvjackclass::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("実行時オプションエラー！");
  nc = id.nc();                 // ケース数変更用
  for (ldata_itr dpos = datalst.begin(); dpos != datalst.end(); ++dpos)
  {
    bool blank = false;
    for (size_t j = 0; blank == false && j < id.ni(); ++j)
    {
      if (dpos->fdata(j) == BLANKf) blank = true;
    }
    if (blank)                  // ケースを削除する
    {
      dpos = datalst.erase(dpos);
      --dpos;  --nc; // 削除した次の要素の位置を返すので1つ戻す必要がある
    }
  }
  if (nc < (int)(id.nc())) id.nc(nc); // ケース数を減らす
  return nc;
}

void cvjackclass::put(option& opt) // 集計結果を出力する
{
  ldata_itr dpos0 = datalst.begin();
  int np = (int)(dpos0->fdata(0)); // パラメータ数
  int nd = (id.ni() - np) / 2; // 予測数
  varclass var(np, nd);        // 変数オブジェクト
  cal_jack(var);
  put_jack(cout, opt, var);
  cout << "\n" << endl;
}

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