// rmat.C Version 7.45  Copyleft (c) 杉本典夫
// Time-stamp: <2026-04-10 13:08:23 tac>
//
// 項目間の相関行列とその検定を計算する。
//
// 使用ライブラリ:libdans.so, libfun.so
// コンパイル法 :g++ -O -s -W rmat rmat.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 =            // プログラム使用法表示用
{
  "rmat  Version 7.45           Copyleft(Vo\\oV) 2026 杉本典夫",
  "書式:rmat  [<オプション> …]  [<データファイル> …]",
  "機能:相関行列・有意確率・有意記号を出力するんよ。",
  "<オプション>",
  "  -f=<ファイル>:オプションを指定したファイルから読み込むんよ。",
  "  -h[=<コマンド名>]:オンラインマニュアルを表示するよ。(-h:-h=man)",
  "  -res=[-][{def|all|i|l|p|r|r2|stat}][,…]:出力する項目を指定するんよ。",
  "    -:-を付けた項目を出力しない(-res=-:全ての項目を出力しない)",
  "    def:オプション省略時と同じ項目を出力する  all:全ての項目を出力する",
  "    i:有意記号行列  l:多重共線性指標  p:有意確率行列  r:相関行列  r2:寄与率行列",
  "    stat:基礎統計量",
  "    これらの記号は連続して複数指定できるんよ。"
};

const string Title = "=== 項目間の単相関 ==="; // 表題

class option                    // オプションクラス
{
private:
  string v_path;                // プログラムのパス
  string v_pname;               // プログラム名
  lstr reslst;                  // 分析結果出力リスト
  mapsb resmap;                 // 分析結果出力の有無マップ
public:
  option(const char* c_pname = nullptr); // コンストラクタ
  ~option()                              // デストラクタ
  {}
  const string& path() const    // プログラムのパスを返す
  { return v_path; }
  const string& pname() const   // プログラム名を返す
  { return v_pname; }
  bool res(const string& code); // 結果出力フラッグを返す
  void setdef()                 // 省略時オプション指定を設定する
  { getoptdef(*this); }
  int check(size_t ni);         // オプションのチェックとセット
  void set(const string& optstr, char devc = OPTdev);
  // オプション指定によってオプションを設定する
  int set_res();                // -resオプションのセット
};

option::option(const char* c_pname) : v_path(), v_pname(), reslst(), resmap()
// コンストラクタ
{
  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::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] == "f")            // オプション指定ファイル
  { err = fgetopt(*this, opt[1], true); }
  else if (opt[0] == "h")       // オンラインヘルプ表示
  { err = help(pname(), USAGE, opt[1]); }
  else if (opt[0] == "res")     // 結果の出力指定
  {
    if (strsplit(reslst, opt[1], 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(), "項目が2つ以上必要だよ！");
  if (set_res() < 0) return -2; // -resオプションのエラー
  return 0;                    // エラーなし
}

int option::set_res()           // -resオプションのセット
{
  resmap["res"] = resmap["def"] = resmap["i"] = resmap["p"] = resmap["r"] = resmap["stat"] = true;
  resmap["r2"] = resmap["l"] = false;
  if (reslst.size() > 0)
  {
    resmap["def"] = resmap["i"] = resmap["p"] = resmap["r"] = 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["i"] = resmap["p"] = resmap["r"] = resmap["stat"] = true; }
      else if (resstr == "all")
      { resmap["def"]= resmap["i"] = resmap["p"] = resmap["r"] = resmap["r2"] = resmap["stat"] = resmap["l"] = true; }
      else if (resitr == resmap.end())
      { return disperr(pname(), "-res=" + resstr + "オプション指定が不適切！"); }
      else
      { resmap[resstr] = true; }
    }
  }
  return 0;                     // エラーなし
}

class varclass                  // 変数クラス
{
private:
  mat<double> v_smat;           // 積和行列
  vecdbl v_mevec;               // 平均値ベクトル
  vecdbl v_vrvec;               // 分散ベクトル
  int v_ni;                     // 項目数
  int v_n;                      // 例数
  restype v_resvar;             // 処理結果の種類
public:
  varclass(size_t ni) : v_smat(ni + 1), v_mevec(ni + 1, 0.0), v_vrvec(ni + 1, 0.0), v_ni(ni), v_n(0), v_resvar(resnon)
  // コンストラクタ
  {}
  ~varclass()                   // デストラクタ
  {}
  const mat<double>& smat() const // 積和行列を返す
  { return v_smat; }
  const vecdbl& mevec() const   // 平均値ベクトルを返す
  { return v_mevec; }
  const vecdbl& vrvec() const   // 分散ベクトルを返す
  { return v_vrvec; }
  int ni() const                // 項目数を返す
  { return v_ni; }
  int n() const                 // 例数を返す
  { return v_n; }
  int df() const                // 自由度を返す
  { return n() - 2; }
  double r(int i, int j) const  // 相関係数を返す
  { return (0 < i && i <= ni() && 0 < j && j <= ni()) ? v_smat(i, j) : 0.0; }
  restype resvar() const        // 処理結果の種類を返す
  { return v_resvar; }
  void cal(option& opt, ldata& datalst); // データを集計して各種統計量を求める
  void s2r();                   // 積和行列を相関行列に変換する
};

void varclass::cal(option& opt, ldata& datalst)
// データを集計して各種統計量を求める
{
  int dim = ni() + 1;
  vecdbl fdatavec(dim, 0.0);
  fdatavec[0] = 1.0;
  for (auto& cdatai : datalst)
  {
    bool blank = false;
    int j = 1;
    while (blank == false && j < dim)
    {
      if (cdatai.fdata(j - 1) == BLANKf)
      { blank = true; }
      else
      {
        fdatavec[j] = cdatai.fdata(j - 1);
        ++j;
      }
    }
    if (blank == false)
    {
      for (int i = 0; i < dim; ++i)
      {
        double x = fdatavec[i];
        for (int j = i; j < dim; ++j) v_smat.add(i, j, x * fdatavec[j]);
      }
    }
  }
  v_n = (int)(v_smat(0, 0));
  for (int i = 0; i < ni(); ++i)
  {
    for (int j = i + 1; j <= ni(); ++j) v_smat.assign(j, i, v_smat(i, j));
  }
  v_smat.sweep(0);
  for (int i = 0; i <= ni(); ++i)
  {
    v_mevec[i] = v_smat(0, i);
    v_vrvec[i] = (v_n > 1) ? v_smat(i, i) / (double)(v_n - 1) : 0.0;
  }
  v_resvar = (v_n > 0) ? resok : reserr;
}

void varclass::s2r()            // 積和行列を相関行列に変換する
{
  smat2rmat(v_smat);
  v_resvar = (n() > 2) ? resok : reserr;
}

class rmatclass                 // rmatメインクラス
{
private:
  string fname;                 // ファイル名
  idd id;                       // データ識別情報
  ldata datalst;                // データリスト
  void put_id(ostream& fout);   // タイトルとデータ識別情報を出力する
  void put_r2mat(ostream& fout, mattype mtype, const varclass& var) const;
  // 寄与率行列を出力する
  void put_pmat(ostream& fout, mattype mtype, const varclass& var) const;
  // 有意確率・有意記号行列を出力する
  void put_ml(ostream& fout, const varclass& var) const;
  // 多重共線性指標を出力する
public:
  rmatclass() : fname(), id(), datalst() // コンストラクタ
  {}
  ~rmatclass()                  // デストラクタ
  {}
  int get(const char* c_fname, option& opt); // データを読み込み累積する
  void put(option& opt);                     // 集計結果を出力する
};

void rmatclass::put_id(ostream& fout) // タイトルとデータ識別情報を出力する
{
  puttitle(fout, Title);
  fout << "\n";
  id.put_title(fout);
  id.put_stra(fout);
  for (size_t i = 0; i < id.ni(); ++i)
  {
    fout << "変数x" << setw(3) << i + 1;
    fout << ":" << id.item(i);
    fout << "\n";
  }
}

void rmatclass::put_r2mat(ostream& fout, mattype mtype, const varclass& var) const
// 寄与率行列を出力する
{
  constexpr int Cmax = 9;       // 行列表横の最大数
  int prec0 = fout.precision(); // 現在の有効桁数を保存する
  int low = 1;
  fout << "\n・寄与率行列(coefficient of determination matrix)";
  while (low <= var.ni())
  {
    int up = (low + Cmax - 1 < var.ni()) ? low + Cmax - 1 : var.ni();
    fout << "\n    ";
    for (int i = low; i <= up; ++i) fout << "    x" << setw(3) << i;
    fout << "\n";
    int llen = 4 + (up - low + 1) * 8;
    putline(fout, llen);
    for (int i = 1; i <= var.ni(); ++i)
    {
      fout << "x" << setw(3) << i;
      fout.precision(4);  fout.setf(ios::fixed);
      for (int j = low; j <= up; ++j)
      {
        if (i == j) fout << "    1   ";
        else fout << setw(8) << var.r(i, j) * var.r(i, j);
      }
      fout.unsetf(ios::fixed);
      fout << "\n";
    }
    fout.precision(prec0);      // 小数点以下の精度を元に戻す
    putline(fout, llen);
    low += Cmax;
  }
}

void rmatclass::put_pmat(ostream& fout, mattype mtype, const varclass& var) const
// 有意確率・有意記号行列を出力する
{
  constexpr int Cmax = 9;       // 行列表横の最大数
  int prec0 = fout.precision(); // 現在の有効桁数を保存する
  int low = 1;
  fout << "\n・有意";
  if (mtype == matr) fout << "確率行列(p値)";
  else fout << "記号行列(・:p>0.1 +:p<0.1 *:p<0.05 **:p<0.01 ***:p<0.001)";
  while (low <= var.ni())
  {
    int up = (low + Cmax - 1 < var.ni()) ? low + Cmax - 1 : var.ni();
    fout << "\n    ";
    for (int i = low; i <= up; ++i) fout << "    x" << setw(3) << i;
    fout << "\n";
    int llen = 4 + (up - low + 1) * 8;
    putline(fout, llen);
    for (int i = 1; i <= var.ni(); ++i)
    {
      fout << "x" << setw(3) << i;
      if (mtype == matr) { fout.precision(4);  fout.setf(ios::fixed); }
      for (int j = low; j <= up; ++j)
      {
        if (i == j)
        { fout << "    -   "; }
        else
        {
          double r = var.r(i, j);
          double f = (sgn(1.0 - r * r) == 0) ? 0.0 : r * r / (1.0 - r * r) * (double)(var.df());
          ldbl pr = ftop(f, 1, var.df());
          if (mtype == matr)
          { fout << setw(8) << pr; }
          else
          {
            if (pr > 0.1L) fout << "    ・  ";
            else fout << "    " << left << setw(4) << pind(pr) << right;
          }
        }
      }
      if (mtype == matr) { fout.precision(prec0);  fout.unsetf(ios::fixed); }
      fout << "\n";
    }
    putline(fout, llen);
    low += Cmax;
  }
}

void rmatclass::put_ml(ostream& fout, const varclass& var) const
// 多重共線性指標を出力する
{
  mat<double> rrmat(var.smat());
  mat<double> tmat(rrmat.sub(1, 1, var.ni(), var.ni()));
  tmat.inv();
  for (int i = 0; i < var.ni(); ++i)
  {
    for (int j = 0; j < var.ni(); ++j) rrmat.assign(1 + i, 1 + j, tmat(i, j));
  }
  putmat(fout, "・相関行列の逆行列", mats, var.n(), rrmat);
  constexpr int LENline = 40;
  fout << "\n・多重共線性の指標";
  fout << "\n変数         VIF      許容度  重相関係数\n";
  putline(fout, LENline);
  for (int i = 0; i < var.ni(); ++i)
  {
    double rii = tmat(i, i);
    double tii = (fabs(rii) > FDeps) ? 1.0 / rii : 0.0;
    fout << "x" << setw(3) << 1 + i << setw(12) << rii << setw(12) << tii;
    fout << setw(12) << sqrt(1.0L - tii) << "\n";
  }
  putline(fout, LENline);
}

int rmatclass::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 rmatclass::put(option& opt) // 集計結果を出力する
{
  if (opt.res("res") == false) return;
  varclass var(id.ni());        // 変数オブジェクト
  var.cal(opt, datalst);        // 各種の計算実行
  put_id(cout);
  if (var.resvar() == reserr)
  {
    cout << "※相関行列計算不可能！\n\n" << endl;
    return;
  }
  if (opt.res("stat")) putmsdvec(cout, var.n(), var.mevec(), var.vrvec());
  if (opt.res("r")) putmat(cout, "・相関行列(correlation coefficient matrix)", matr, var.n(), var.smat());
  var.s2r();                    // 積和行列を相関行列に変換する
  if (var.resvar() != reserr)
  {
    if (opt.res("r2")) put_r2mat(cout, matetc, var);
    if (opt.res("p")) put_pmat(cout, matr, var);
    if (opt.res("i")) put_pmat(cout, matetc, var);
    if (opt.res("l")) put_ml(cout, var);
  }
  cout << "\n" << endl;
}

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