// cca.C Version 7.44  Copyleft (c) 杉本典夫
// Time-stamp: <2025-10-30 17:04:50 tac>
//
// 指定した2つの項目グループ間で正準相関分析を行う。
//
// 使用ライブラリ:libdans.so, libfun.so
// コンパイル法 :g++ -O -s -W cca cca.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 =            // プログラム使用法表示用
{
  "cca  Version 7.44            Copyleft(Vo\\oV) 2025 杉本典夫",
  "書式:cca [<オプション> …]  [<データファイル> …]",
  "機能:指定した2つの項目グループ間で正準相関分析をするのだ。",
  "<オプション>",
  "  -dnsname[=<出力ファイル名>]:DANSデータを出力するファイル名を指定するのだ。",
  "    オプション省略時:DANSデータ出力せず  -dnsname:-dnsname=cca_dns.d",
  "    出力項目:正準スコア",
  "  -f=<ファイル>:指定したファイルからオプションを読み込むのだ。",
  "  -h[=<コマンド名>]:オンラインマニュアルを表示するのだ。(-h:-h=man)",
  "  -res=[-][{def|all|ccc|ccl|cl|cvc|stat}][,…]:結果の出力項目を指定するのだ。",
  "    -:-を付けた項目を出力しない(-res=-:全ての項目を出力しない)",
  "    def:オプション省略時と同じ項目を出力する  all:全ての項目を出力する",
  "    ccc:正準相関係数  ccl:交差負荷量  cl:正準負荷量  cvc:正準変量係数",
  "    stat:基礎統計量  オプション省略時:全て出力",
  "  -x=[<開始項目番号>][-[<終了項目番号>]]:変数Xの項目番号を指定するのだ。",
  "    オプション省略時:-x=1",
  "  -y=[<開始項目番号>][-[<終了項目番号>]]:変数Yの項目番号を指定するのだ。",
  "    オプション省略時:-y=2"
};

const string Title = "=== 正準相関分析(canonical correlation analysis) ==="; // 表題
const string DNSNAMEdef = "cca_dns.d"; // 省略時DANSファイル名

class option                    // オプションクラス
{
private:
  string v_path;                // プログラムのパス
  string v_pname;               // プログラム名
  bool v_dns;                   // 正準スコアの出力フラッグ
  string v_dnsname;             // DANSファイル名
  lstr reslst;                  // 分析結果出力リスト
  mapsb resmap;                 // 分析結果出力の有無マップ
  bool v_stat;                  // 基礎統計量の出力フラッグ
  pairint xno;                  // 変数xの開始・終了項目番号
  pairint yno;                  // 変数yの開始・終了項目番号
public:
  option(const char* c_pname = nullptr); // コンストラクタ
  ~option()                              // デストラクタ
  {}
  const string& path() const    // プログラムのパスを返す
  { return v_path; }
  const string& pname() const   // プログラム名を返す
  { return v_pname; }
  bool dns() const              // 正準スコアの出力フラッグを返す
  { return v_dns; }
  const string& dnsname() const // DANSデータファイル名を返す
  { return v_dnsname; }
  bool stat() const             // 基礎統計量出力フラッグを返す
  { return v_stat; }
  int xs() const                // 変数xの開始番号を返す
  { return xno.first; }
  int xe() const                // 変数xの終了番号を返す
  { return xno.second; }
  int nx() const                // 変数xの数を返す
  { return xno.second - xno.first + 1; }
  int ys() const                // 変数yの開始番号を返す
  { return yno.first; }
  int ye() const                // 変数yの終了番号を返す
  { return yno.second; }
  int ny() const                // 変数yの数を返す
  { return yno.second - yno.first + 1; }
  bool res(const string& code); // 結果出力フラッグを返す
  void setdef()                 // 省略時オプション指定を設定する
  { getoptdef(*this); }
  void set(const string& optstr, char devc = OPTdev);
  // オプション指定によってオプションを設定する
  int check(size_t ni);         // オプションのチェックとセット
  int set_res();                // -resオプションのセット
};

option::option(const char* c_pname) : v_path(), v_pname(), v_dns(false), v_dnsname(), reslst(), resmap(), v_stat(true), xno(1, 1), yno(2, 2)
// コンストラクタ
{
  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] == "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] == "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 if (opt[0] == "x")       // 変数Xの項目番号指定
  {
    xno.first = 1;  xno.second = INT_MAX;
    if (getrange(xno.first, xno.second, opt[1]) < 1) err = badopt(pname(), optstr, USAGE);
    else if (xno.first < 1 || xno.first > xno.second) err = badopt(pname(), optstr, USAGE);
  }
  else if (opt[0] == "y")       // 変数yの項目番号指定
  {
    yno.first = 1;  yno.second = INT_MAX;
    if (getrange(yno.first, yno.second, opt[1]) < 1) err = badopt(pname(), optstr, USAGE);
    else if (yno.first < 1 || yno.first > yno.second) 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 ((int)ni < xno.second || (int)ni < yno.second) return disperr(pname(), "項目番号指定が不適なのだ！");
  if (set_res() < 0) return -3; // -resオプションのエラー
  return 0;                     // エラーなし
}

int option::set_res()           // -resオプションのセット
{
  resmap["res"] = resmap["def"] = resmap["ccc"] = resmap["ccl"] = resmap["cl"] = resmap["cvc"] = resmap["stat"] = true;
  if (reslst.size() > 0)
  {
    resmap["def"] = resmap["ccc"] = resmap["ccl"] = resmap["cl"] = resmap["cvc"] = 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" || resstr == "all") // 全ての結果を出力する
      {
        resmap["def"] = resmap["ccc"] = resmap["ccl"] = resmap["cl"] = resmap["cvc"] = resmap["stat"] = true;
      }
      else if (resitr == resmap.end())
      { return disperr(pname(), "-res=" + resstr + "オプション指定が不適切！"); }
      else
      { resmap[resstr] = true; }
    }
  }
  return 0;                     // エラーなし
}

class varclass                 // 変数クラス
{
private:
  mat<ldbl> v_smat;             // 積和行列
  int v_n;                      // 例数
  int v_p;                      // パラメータ数
  int v_q;                      // 因子数
  vecldbl v_mevec;              // 平均値ベクトル
  vecldbl v_vrvec;              // 分散ベクトル
  vecldbl evalvec;              // 固有値ベクトル
  mat<ldbl> v_xcvcmat;          // Xの正準変量係数行列
  mat<ldbl> v_xclmat;           // Xの正準負荷量行列
  mat<ldbl> v_xcclmat;          // Xの交差負荷量行列
  mat<ldbl> v_ycvcmat;          // Yの正準変量係数行列
  mat<ldbl> v_yclmat;           // Yの正準負荷量行列
  mat<ldbl> v_ycclmat;          // Yの交差負荷量行列
  int cal_smat(option& opt, ldata& datalst); // 積和行列を求める
  int cal_cca(option& opt);    // 正準相関を計算する
  void evmat_cvcmat(mat<ldbl>& cvcmat, int is, int ie);
  // 固有ベクトルを正準変量係数に変換する
  void evmat_clmat(mat<ldbl>& clmat, int is, int ie);
  // 固有ベクトルを正準負荷量(正準変量との相関係数)に変換する
public:
  varclass(option& opt);        // コンストラクタ
  ~varclass()                   // デストラクタ
  {}
  int n() const                 // 例数を返す
  { return v_n; }
  int p() const                 // パラメータ数を返す
  { return v_p; }
  int q() const                 // 因子数を返す
  { return v_q; }
  const mat<ldbl>& smat() const // 積和行列を返す
  { return v_smat; }
  const vecldbl& mevec() const  // 平均値ベクトルを返す
  { return v_mevec; }
  const vecldbl& vrvec() const  // 分散ベクトルを返す
  { return v_vrvec; }
  ldbl eval(int idx) const      // 固有値を返す
  { return (0 <= idx && idx <= q()) ? evalvec[idx] : 0.0L; }
  const mat<ldbl>& xcvcmat() const // Xの正準変量係数行列を返す
  { return v_xcvcmat; }
  mat<ldbl>& xclmat()           // Xの正準負荷量行列を返す
  { return v_xclmat; }
  mat<ldbl>& xcclmat()          // Xの交差負荷量行列を返す
  { return v_xcclmat; }
  const mat<ldbl>& ycvcmat() const // Yの正準変量係数行列を返す
  { return v_ycvcmat; }
  mat<ldbl>& yclmat()           // Yの正準負荷量行列を返す
  { return v_yclmat; }
  mat<ldbl>& ycclmat()          // Yの交差負荷量行列を返す
  { return v_ycclmat; }
  void cal(option& opt, ldata& datalst); // データを集計して各種統計量を求める
  void calscore(option& opt, datac& score, datac& cdata) const;
  // 正準スコアを計算して返す
};

varclass::varclass(option& opt) : v_smat(), v_n(0), v_p(0), v_q(0), v_mevec(), v_vrvec(), evalvec(), v_xcvcmat(), v_xclmat(), v_xcclmat(), v_ycvcmat(), v_yclmat(), v_ycclmat()
// コンストラクタ
{
  v_p = opt.nx() + opt.ny();    // パラメータ数
  size_t dim = 1 + v_p;         // 行列とベクトルの次元数
  v_smat.resize(dim, dim);
  v_mevec.resize(dim, 0.0);  v_vrvec.resize(dim, 0.0);
  int qmax = min(opt.nx(), opt.ny());
  evalvec.resize(qmax + 1, 0.0);
  v_xcvcmat.resize(opt.nx() + 1, qmax + 1);
  v_xclmat.resize(opt.nx() + 1, qmax + 1);
  v_xcclmat.resize(opt.nx() + 1, qmax + 1);
  v_ycvcmat.resize(opt.ny() + 1, qmax + 1);
  v_yclmat.resize(opt.ny() + 1, qmax + 1);
  v_ycclmat.resize(opt.ny() + 1, qmax + 1);
}

int varclass::cal_smat(option& opt, ldata& datalst)
// 積和行列を求める
{
  size_t dim = 1 + p();
  vecdbl fdatavec(dim);
  v_smat = 0.0L;  fdatavec[0] = 1.0L;
  for (auto& cdatai : datalst)
  {
    bool blank = false;
    for (int j = opt.xs(); blank == false && j <= opt.xe(); ++j)
    {
      int k = 1 + j - opt.xs();
      fdatavec[k] = cdatai.fdata(j - 1);
      if (fdatavec[k] == BLANKf) blank = true;
    }
    for (int j = opt.ys(); blank == false && j <= opt.ye(); ++j)
    {
      int k = 1 + opt.nx() + j - opt.ys();
      fdatavec[k] = cdatai.fdata(j - 1);
      if (fdatavec[k] == BLANKf) blank = true;
    }
    if (blank == false)
    {
      for (size_t i = 0; i < dim; ++i)
      {
        double x = fdatavec[i];
        for (size_t j = i; j < dim; ++j) v_smat.add(i, j, (ldbl)(x * fdatavec[j]));
      }
    }
  }
  for (size_t i = 0; i < dim - 1; ++i)
  {
    for (size_t j = i + 1; j < dim; ++j) v_smat.assign(j, i, v_smat(i, j));
  }
  v_n = (int)(v_smat(0, 0));
  v_smat.sweep(0);              // 単純積和行列を積和行列に変換する
  return v_n;                   // 例数を返す
}

int varclass::cal_cca(option& opt)
// 正準相関を計算する
{
  v_q = 0;
  if (v_n < p()) return v_q;    // 正順相関係数計算不可能
  mat<ldbl> rmat = v_smat.sub(1, 1, v_smat.r() - 1, v_smat.c() - 1);
  smat2rmat(rmat);
  mat<ldbl> rxxmat = rmat.sub(0, 0, opt.nx() - 1, opt.nx() - 1);
  mat<ldbl> rxxrmat = rxxmat;
  rxxrmat.inv();
  mat<ldbl> rxymat = rmat.sub(0, opt.nx(), opt.nx() - 1, p() - 1);
  mat<ldbl> ryyrmat = rmat.sub(opt.nx(), opt.nx(), p() - 1, p() - 1);
  ryyrmat.inv();
  mat<ldbl> ryxmat = rmat.sub(opt.nx(), 0, p() - 1, opt.nx() - 1);
  mat<ldbl> tmat = rxxrmat * rxymat * ryyrmat * ryxmat;
  evalvec[0] = 0.0;
  for (int i = 0; i < opt.nx(); i++) evalvec[0] += tmat(i, i);
  vec<ldbl> axv(opt.nx()), ayv(opt.nx());
  int qmax = min(opt.nx(), opt.ny());
  while (v_q < qmax && (evalvec[v_q + 1] = eigenud(tmat, axv, ayv)) > 0.0L)
  {
    ++v_q;
    vec<ldbl> tv = rxxmat * axv;
    ldbl w = inp(axv, tv);
    axv *= (w > 0.0L) ? sqrt(1.0L / w) : 0.0L;
    vec<ldbl> bv = ryyrmat * ryxmat * axv;
    bv *= sqrt(1.0L / evalvec[v_q]);
    for (int i = 0; i < opt.nx(); ++i) v_xcvcmat.assign(i + 1, v_q, axv(i));
    for (int i = 0; i < opt.ny(); ++i) v_ycvcmat.assign(i + 1, v_q, bv(i));
  }
  v_xcvcmat.resize(opt.nx() + 1, v_q + 1);
  v_ycvcmat.resize(opt.ny() + 1, v_q + 1);
  return v_q;                   // 正順相関係数の数を返す
}

void varclass::evmat_cvcmat(mat<ldbl>& cvcmat, int is, int ie)
// 固有ベクトルを正準変量係数に変換する
{
  int p = ie - is + 1, q = cvcmat.c() - 1;
  vecldbl sdrvec(p, 0.0L);
  for (int i = is; i <= ie; ++i)
  { sdrvec[i - is] = (v_smat(i, i) > 0.0L) ? sqrt((ldbl)(n() - 1) / v_smat(i, i)) : 0.0L; }
  for (int j = 0; j <= q; ++j) cvcmat.assign(0, j, 0.0L);
  for (int i = 1; i <= p; ++i)
  {
    ldbl w = sdrvec[i - 1], m = v_smat(0, i + is - 1);
    for (int j = 1; j <= q; ++j)
    {
      cvcmat.assign(i, j, cvcmat(i, j) * w);
      cvcmat.add(0, j, -cvcmat(i, j) * m);
    }
  }
}

void varclass::evmat_clmat(mat<ldbl>& clmat, int is, int ie)
// 固有ベクトルを正準負荷量(正準変量との相関係数)に変換する
{
  mat<ldbl> rmat = v_smat.sub(is, is, ie, ie);
  smat2rmat(rmat);
  mat<ldbl> tmat = clmat.sub(1, 1, clmat.r() - 1, clmat.c() - 1);
  mat<ldbl> lmat = rmat * tmat;
  for (size_t i = 1; i < clmat.r(); ++i)
  {
    for (size_t j = 1; j < clmat.c(); ++j) clmat.assign(i, j, lmat(i - 1, j - 1));
  }
}

void varclass::cal(option& opt, ldata& datalst)
// データを集計して各種統計量を求める
{
  if (cal_smat(opt, datalst) < 1) return;
  size_t dim = 1 + p();            // 行列とベクトルの次元数
  for (size_t i = 0; i < dim; ++i)   // 平均値と分散をセットする
  {
    v_mevec[i] = v_smat(0, i);
    v_vrvec[i] = (n() > 1) ? v_smat(i, i) / (ldbl)(n() - 1) : 0.0;
  }
  if (cal_cca(opt) < 1) return;
  v_xclmat = v_xcvcmat;  v_yclmat = v_ycvcmat;
  evmat_cvcmat(v_xcvcmat, 1, opt.nx());
  evmat_cvcmat(v_ycvcmat, opt.nx() + 1, p());
  evmat_clmat(v_xclmat, 1, opt.nx());
  v_xcclmat = v_xclmat;
  for (size_t i = 1; i < v_xcclmat.r(); ++i)
  {
    for (size_t j = 1; j < v_xcclmat.c(); ++j) v_xcclmat.assign(i, j, v_xcclmat(i, j) * sqrt(evalvec[j]));
  }
  evmat_clmat(v_yclmat, opt.nx() + 1, p());
  v_ycclmat = v_yclmat;
  for (size_t i = 1; i < v_ycclmat.r(); ++i)
  {
    for (size_t j = 1; j < v_ycclmat.c(); ++j) v_ycclmat.assign(i, j, v_ycclmat(i, j) * sqrt(evalvec[j]));
  }
}

void varclass::calscore(option& opt, datac& score, datac& cdata) const
// 正準スコアを計算して返す
{
  vecldbl fdatavec(opt.nx() + opt.ny(), 0.0L);
  bool blank = false;
  for (int j = opt.xs(); blank == false && j <= opt.xe(); ++j)
  {
    int k = j - opt.xs();
    fdatavec[k] = cdata.fdata(j - 1);
    if (fdatavec[k] == BLANKf) blank = true;
  }
  for (int j = opt.ys(); blank == false && j <= opt.ye(); ++j)
  {
    int k = opt.nx() + j - opt.ys();
    fdatavec[k] = cdata.fdata(j - 1);
    if (fdatavec[k] == BLANKf) blank = true;
  }
  if (blank)                    // スコアの計算ができない
  { score = BLANKs; }
  else
  {
    for (int i = 1; i <= q(); ++i)
    {
      ldbl z = v_xcvcmat(0, i);
      for (int j = 0; j < opt.nx(); ++j) z += v_xcvcmat(j + 1, i) * fdatavec[j];
      score.assign(i - 1, z);
    }
    for (int i = 1; i <= q(); ++i)
    {
      ldbl z = v_ycvcmat(0, i);
      for (int j = 0; j < opt.ny(); ++j) z += v_ycvcmat(j + 1, i) * fdatavec[opt.nx() + j];
      score.assign(q() + i - 1, z);
    }
  }
}

class ccaclass                  // ccaメインクラス
{
private:
  string fname;                 // ファイル名
  idd id;                       // データ識別情報
  ldata datalst;                // データリスト
  void put_resid(ostream& fout, option& opt);
  // 解析結果のタイトルとデータ識別情報を出力する
  void put_ccatab(ostream& fout, option& opt, varclass& var);
  // 正準相関分析結果を出力する
  void put_cl(ostream& fout, const string& mtitle, const string& rhed, const string& ched, mat<ldbl>& lmat) const;
  // 正準負荷量または交差負荷量を出力する。
  void put_dnscca(option& opt, const varclass& var);
  // 正準スコアをデータファイル形式で出力する
  void cal_score(option& opt, datac& score, datac& cdata, const mat<ldbl>& amat, const mat<ldbl>& bmat, int q) const;
  // 正準スコアを計算し、データとして返す
public:
  ccaclass() : fname(), id(), datalst() // コンストラクタ
  {}
  ~ccaclass()                   // デストラクタ
  {}
  int get(const char* c_fname, option& opt); // データを読み込み累積する
  void put(option& opt);                     // 集計結果を出力する
};

void ccaclass::put_resid(ostream& fout, option& opt)
// 解析結果のタイトルとデータ識別情報を出力する
{
  puttitle(fout, Title);
  fout << "\n";
  id.put_title(fout);
  id.put_stra(fout);
  for (int i = opt.xs(); i <= opt.xe(); ++i)
  {
    fout<< "変数x" << setw(3) << i - opt.xs() + 1;
    fout << ":" << id.item(i - 1) << "\n";
  }
  for (int i = opt.ys(); i <= opt.ye(); ++i)
  {
    fout << "変数y" << setw(3) << i - opt.ys() + 1;
    fout << ":" << id.item(i - 1) << "\n";
  }
}

void ccaclass::put_ccatab(ostream& fout, option& opt, varclass& var)
// 正準相関分析結果を出力する
{
  constexpr int LENline = 72;   // 区切り線の長さ
  int p = opt.nx() + opt.ny();
  if (opt.res("ccc"))
  {
    fout << "\n・正準相関分析結果\n";
    fout << "番号    正準相関係数  固有値(寄与率)     正規分布z値     有意確率p値\n";
    putline(fout, LENline);
    ldbl plamda = 1.0;
    for (int i = 1; i <= var.q(); ++i)
    {
      ldbl lamda = var.eval(i), z = 0.0L;
      if (lamda >= 1.0)
      { lamda = 1.0; }
      else
      {
        plamda *= (1.0L - lamda);
        z = sqrt((ldbl)(var.n()) * lamda / (1.0L - lamda));
      }
      ldbl pr = ztop(z);
      fout << "No" << setw(2) << i << setw(16) << sqrt(lamda);
      fout << setw(16) << lamda << setw(16) << z;
      fout << setw(16) << pr << pind(pr) << "\n";
    }
    putline(fout, LENline);
    ldbl csqr = - ((ldbl)(var.n() - 1) - (ldbl)(p + 1) * 0.5L) * log(plamda);
    int df = opt.nx() * opt.ny();
    ldbl pr = ctop(csqr, df);
    fout << "全ての正準相関係数=0の検定:χ2乗=" << csqr << "  自由度=" << df << "  有意確率p=" << pr << pind(pr) << "\n";
  }
  if (opt.res("cvc"))
  {
    putcoef(fout, "・正準変量係数", "U", "x", var.xcvcmat());
    putcoef(fout, "・正準変量係数", "V", "y", var.ycvcmat());
  }
  if (opt.res("cl")) put_cl(fout, "・正準負荷量(変数xと正準変量Uの相関係数)", "U", "x", var.xclmat());
  if (opt.res("ccl")) put_cl(fout, "・交差負荷量(変数xと正準変量Vの相関係数)", "V", "x", var.xcclmat());
  if (opt.res("cl")) put_cl(fout, "・正準負荷量(変数yと正準変量Vの相関係数)", "V", "y", var.yclmat());
  if (opt.res("ccl")) put_cl(fout, "・交差負荷量(変数yと正準変量Uの相関係数)", "U", "y", var.ycclmat());
}

void ccaclass::put_cl(ostream& fout, const string& mtitle, const string& rhed, const string& ched, mat<ldbl>& lmat) const
// 正準負荷量または交差負荷を量出力する。
{
  constexpr size_t Cmax = 6;    // 行列表横の最大数
  size_t pidx = lmat.r() - 1, qidx = lmat.c() - 1, low = 1, up = 0;
  vecldbl l2vec(lmat.c(), 0.0L); // 負荷量平方和
  bool eop = false;             // 表出力終了フラッグ
  for (size_t i = 0; i <= pidx; ++i) lmat.assign(i, 0, 0.0L);
  fout << "\n" << mtitle;
  while (eop == false)
  {
    size_t nc = low + Cmax - 1;
    if (low > qidx) eop = true;
    else if (nc > qidx) { up = qidx;  eop = true; }
    else if (nc == qidx) up = qidx;
    else up = nc;
    fout << "\n変数";
    for (size_t j = low; j <= up; ++j)
    {
      int npad = rhed.size() - strclen(rhed, CHARwide);
      fout << setw(9 + npad) << rhed << setw(3) << j;
    }
    if (eop) fout << "      寄与率";
    fout << "\n";
    size_t llen = 4 + 12 * (up - low + 1 + (eop ? 1 : 0));
    putline(fout, llen);
    for (size_t i = 1; i <= pidx; ++i)
    {
      fout << setw(1) << ched << setw(3) << i;
      int prec0 = fout.precision();
      fout.precision(6);
      fout.setf(ios::fixed);
      for (size_t j = low; j <= up; ++j)
      {
        fout << setw(12) << lmat(i, j);
        ldbl l2 = lmat(i, j) * lmat(i, j);
        lmat.assign(i, 0, lmat(i, 0) + l2);
        l2vec[j] += l2;
      }
      if (eop) fout << setw(12) << lmat(i, 0);
      fout << "\n";
      fout.unsetf(ios::fixed);
      fout.precision(prec0);
    }
    putline(fout, llen);
    int prec0 = fout.precision();
    fout.precision(6);
    fout.setf(ios::fixed);
    fout << "寄与";
    for (size_t j = low; j <= up; ++j)
    {
      ldbl r2 = l2vec[j] / (ldbl)pidx;
      fout << setw(12) << r2;
      l2vec[0] += r2;
    }
    if (eop) fout << setw(12) << l2vec[0];
    fout << "\n";
    low = up + 1;
  }
}

void ccaclass::put_dnscca(option& opt, const varclass& var)
// 正準スコアをデータファイル形式で出力する
{
  ofstream fout;
  if (ofopen(fout, (opt.dnsname()).c_str()) == false) return; // ファイルオープン失敗
  id.put_title(fout, formtag);
  id.put_stra(fout, formtag);
  idd tid(id.title(), var.q() * 2);
  ostringstream ostr;
  for (int i = 0; i < var.q(); ++i)
  {
    ostr << "正準スコアU" << i + 1;
    tid.assign_item(i, ostr.str());
    ostr.str("");
  }
  for (int i = 0; i < var.q(); ++i)
  {
    ostr << "正準スコアV" << i + 1;
    tid.assign_item(var.q() + i, ostr.str());
    ostr.str("");
  }
  tid.put_item(fout, formtag);
  datac score(var.q() * 2);
  for (auto& cdatai : datalst)
  {
    var.calscore(opt, score, cdatai);
    score.put(fout);
  }
}

int ccaclass::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 ccaclass::put(option& opt) // 集計結果を出力する
{
  varclass var(opt);            // 変数オブジェクト
  var.cal(opt, datalst);        // 各種の計算実行
  if (opt.res("res"))
  {
    put_resid(cout, opt);
    if (var.n() < 1)
    { cout << "※有効なデータなし！\n"; }
    else
    {
      if (opt.res("stat"))
      {
        putmsdvec(cout, var.n(), var.mevec(), var.vrvec(), opt.nx(), opt.ny());
        putmat(cout, "・相関行列", matr, var.n(), var.smat(), opt.nx(), opt.ny());
      }
      if (var.q() < 1) cout << "※正順相関係数なし！\n";
      else put_ccatab(cout, opt, var);
    }
    cout << "\n" << endl;
  }
  if (var.q() > 0 && opt.dnsname() != "") put_dnscca(opt, var);
}

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