// fact.C Version 7.44  Copyleft (c) 杉本典夫
// Time-stamp: <2025-10-31 08:57:51 tac>
//
// 指定した項目について因子分析を行う。
//
// 使用ライブラリ:libdans.so, libfun.so
// コンパイル法 :g++ -O -s -W fact fact.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 =            // プログラム使用法表示用
{
  "fact  Version 7.44           Copyleft(Vo\\oV) 2025 杉本典夫",
  "書式:fact  [<オプション> …]  [<データファイル> …]",
  "機能:全項目を対象とした因子分析をするよ。",
  "<オプション>",
  "  -f=<ファイル>:指定したファイルからオプションを読み込むよ。",
  "  -d2[={max|q|<係数>}]:独自性の初期値を指定するよ。(0＜係数＜1)",
  "  -dnsname[=<出力ファイル名>]:DANSデータを出力するファイル名を指定します。",
  "    オプション省略時:DANSデータ出力せず  -dnsname:-dnsname=fact_dns.d",
  "    出力項目:因子スコア",
  "  -eps=<収束誤差>:収束誤差を指定するよ。(0＜収束誤差)",
  "  -h[=<コマンド名>]:オンラインマニュアルを表示するよ。(-h:-h=man)",
  "  -lamda=<基準固有値>:因子抽出時の基準固有値を指定するよ。(0≦基準固有値)",
  "  -meth={ml|pf}:因子抽出法の指定だよ。",
  "    ml:最尤法  pf:主因子法(オプション省略時)",
  "  -ngpname[=<出力ファイル名>]:Ngraph用データの出力ファイル名を指定します。",
  "    オプション省略時:Ngraph用データ出力せず  -ngpname:-ngpname=fact_ngp.dat",
  "    出力項目:因子負荷量",
  "  -pm[=<ベキ乗数>]:promax回転の指定だよ。(-pm:pm=4)",
  "  -q=<因子数>:因子抽出時の最大因子数を指定するよ。(0＜因子数≦変数の数)",
  "  -res=[-][{def|all|fa|fc|fl|fr|ft|lt|stat}][,…]",
  "    :解析結果の出力項目を指定するよ。",
  "    -:-を付けた項目を出力しない(-res=-:全ての項目を出力しない)",
  "    def:オプション省略時と同じ項目を出力する  all:全ての項目を出力する",
  "    fa:因子  fc:因子スコア計算用係数  fl:因子負荷量  fr:因子スコア間相関係数",
  "    ft:因子軸回転用変換行列  lt:尤度比検定の結果  stat:基礎統計量",
  "  -rmax=<最大回数>:反復計算の最大回数を指定するよ。(0≦最大回数)",
  "  -sort[={-|+}]:因子負荷量を因子順に並べ替える指定だよ。(-sort:-sort=+)",
  "  -vm[={-|+}]:規準varimax回転の指定だよ。(-vm:-vm=+)"
};

const string Title = "=== 因子分析(factor analysis) ==="; // 表題
const string DNSNAMEdef = "fact_dns.d"; // 省略時DANSファイル名
const string NGPNAMEdef = "fact_ngp.dat"; // 省略時Ngraphファイル名
constexpr ldbl EPSdef = 1.0e-5L;          // 収束誤差省略時値
constexpr int RMAXdef = 200;              // 反復計算の省略時最大回数
constexpr int PMdef = 4;        // プロマックス回転の省略時次数

enum d2type                     // 独自性初期値の種類
{
  d2max,                        // 上限値
  d2q,                          // 最大因子数に基づく推定値
  d2coef,                       // 上限値にかける係数指定
  d2etc                         // その他
};

enum methtype                   // 因子抽出法の種類
{
  methnon,                      // 指定なし
  methml,                       // 最尤法
  methpf,                       // 主因子法
  methetc                       // その他
};

enum rottype                    // 因子軸回転の種類
{
  rotnon,                       // 回転なし
  rotpm,                        // promax回転
  rotvm,                        // 規準Varimax回転
  rotetc                        // その他
};

class option                    // オプションクラス
{
private:
  string v_path;                // プログラムのパス
  string v_pname;               // プログラム名
  d2type v_d2;                  // 独自性初期値の種類
  ldbl v_d2c;                   // 独自性初期値の係数
  string v_dnsname;             // DANSファイル名
  ldbl v_eps;                   // 収束誤差
  ldbl v_lamda;                 // 因子抽出用因子負荷量
  methtype v_meth;              // 因子抽出法
  string v_ngpname;             // Ngraphファイル名
  int v_pm;                     // promax回転のベキ乗数
  int v_q;                      // 因子抽出用因子数
  lstr reslst;                  // 解析結果出力リスト
  mapsb resmap;                 // 解析結果出力の有無マップ
  int v_rmax;                   // 反復計算の最大回数
  bool v_sort;                  // 因子負荷量並べ替えフラッグ
  bool v_vm;                    // 規準varimax回転実施フラッグ
public:
  option(const char* c_pname = nullptr); // コンストラクタ
  ~option()                              // デストラクタ
  {}
  const string& path() const    // プログラムのパスを返す
  { return v_path; }
  const string& pname() const   // プログラム名を返す
  { return v_pname; }
  d2type d2() const             // 独自性初期値の種類を返す
  { return v_d2; }
  ldbl d2c() const              // 独自性初期値の係数を返す
  { return v_d2c; }
  const string& dnsname() const // DANSデータファイル名を返す
  { return v_dnsname; }
  ldbl eps() const              // 収束誤差を返す
  { return v_eps; }
  ldbl lamda() const            // 因子抽出用因子負荷量を返す
  { return v_lamda; }
  methtype meth() const         // 因子抽出法を返す
  { return v_meth; }
  const string& ngpname() const // Ngraphデータファイル名を返す
  { return v_ngpname; }
  int pm() const                // promax回転のベキ乗数を返す
  { return v_pm; }
  int q() const                 // 因子抽出用因子数を返す
  { return v_q; }
  bool res(const string& code); // 分析結果出力フラッグを返す
  int rmax() const              // 反復計算の最大回数を返す
  { return v_rmax; }
  bool sort() const             // 因子負荷量並べ替えフラッグを返す
  { return v_sort; }
  bool vm() const               // 規準varimax回転実施フラッグを返す
  { return v_vm; }
  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_d2(d2max), v_d2c(1.0L), v_dnsname(), v_eps(EPSdef), v_lamda(1.0L), v_meth(methpf), v_pm(0), v_q(0), v_ngpname(), reslst(), resmap(), v_rmax(RMAXdef), v_sort(false), v_vm(true)
// コンストラクタ
{
  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] == "d2")      // 独自性初期値の種類指定
  {
    if (opt[1].empty() || opt[1] == "max")
    { v_d2 = d2max; }
    else if (opt[1] == "q")
    { v_d2 = d2q; }
    else
    {
      v_d2c = strtodf(opt[1]);
      if (v_d2c <= 0.0L || 1.0L <= v_d2c) err = badopt(pname(), optstr, USAGE);
      else v_d2 = d2coef;
    }
  }
  else 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] == "eps")     // 収束誤差の指定
  {
    v_eps = opt[1].empty() ? EPSdef : strtodf(opt[1]);
    if (v_eps < FDeps) err = badopt(pname(), optstr, USAGE, FDeps, FDmax);
  }
  else if (opt[0] == "h")       // オンラインヘルプ表示
  { err = help(pname(), USAGE, opt[1]); }
  else if (opt[0] == "lamda")   // 因子抽出用因子負荷量の指定
  {
    v_lamda = opt[1].empty() ? 1.0L : strtodf(opt[1]);
    if (v_lamda < 0.0L || v_lamda == (ldbl)BLANKf) err = badopt(pname(), optstr, USAGE);
  }
  else if (opt[0] == "meth")    // 因子抽出法指定
  {
    if (opt[1].empty() || opt[1] == "pf") v_meth = methpf;
    else if (opt[1] == "ml") v_meth = methml;
    else err = badopt(pname(), optstr, USAGE);
  }
  else if (opt[0] == "ngpname") // Ngraphファイル名指定
  {
    if (opt[1].empty()) v_ngpname = NGPNAMEdef;
    else v_ngpname = opt[1];
    if (v_ngpname.size() < 3) err = badopt(pname(), optstr, USAGE);
  }
  else if (opt[0] == "pm")      // promax回数のベキ乗数指定
  {
    char* ptr = nullptr;
    v_pm = opt[1].empty() ? PMdef : strtol(opt[1].c_str(), &ptr, 0);
    if (v_pm < 0) err = badopt(pname(), optstr, USAGE, 0.0, 10.0);
  }
  else if (opt[0] == "q")       // 最大因子数指定
  {
    char* ptr = nullptr;
    v_q = strtol(opt[1].c_str(), &ptr, 0);
    if (v_q < 1) err = badopt(pname(), optstr, USAGE, 1.0, (double)INT_MAX);
  }
  else if (opt[0] == "res")     // 分析結果の出力指定
  {
    if (strsplit(reslst, opt[1], CHRdelm) < 1) err = badopt(pname(), optstr, USAGE);
  }
  else if (opt[0] == "rmax")    // 最大反復回数指定
  {
    char* ptr = nullptr;
    v_rmax = strtol(opt[1].c_str(), &ptr, 0);
    if (v_rmax < 0) err = badopt(pname(), optstr, USAGE, 0.0, (double)INT_MAX);
  }
  else if (opt[0] == "sort")    // 因子負荷量並べ替え指定
  {
    if (opt[1].empty() || opt[1] == "+") v_sort = true;
    else if (opt[1] == "-") v_sort = false;
    else err = badopt(pname(), optstr, USAGE);
  }
  else if (opt[0] == "vm")      // 規準varimax回転実施指定
  {
    if (opt[1].empty() || opt[1] == "+") v_vm = true;
    else if (opt[1] == "-") v_vm = false;
    else 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(), "項目数が足らんよ！");
  if (set_res() < 0) return -3; // -resオプションのエラー
  return 0;                     // エラーなし
}

int option::set_res()           // -resオプションのセット
{
  resmap["res"] = resmap["def"] = resmap["fa"] = resmap["fl"] = resmap["fc"] = resmap["fr"] = resmap["lt"] = resmap["stat"] = true;
  resmap["ft"] = false;
  if (reslst.size() > 0)
  {
    resmap["def"] = resmap["fa"] = resmap["fl"] = resmap["fc"] = resmap["fr"] = resmap["ft"] = resmap["lt"] = 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["fa"] = resmap["fl"] = resmap["fc"] = resmap["fr"] = resmap["lt"] = resmap["stat"] = true;
      }
      else if (resstr == "all") // 全て出力
      {
        resmap["def"] = resmap["fa"] = resmap["fl"] = resmap["fc"] = resmap["fr"] = resmap["ft"] = resmap["lt"] = 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;             // 積和行列
  vecldbl v_mevec;              // 平均値ベクトル
  vecldbl v_vrvec;              // 分散ベクトル
  vecldbl evalvec0;             // 初期固有値ベクトル
  vecldbl evalvec;              // 固有値ベクトル
  mat<ldbl> v_lmat;             // 因子負荷量行列
  mat<ldbl> v_amat;             // 因子スコア計算用係数行列
  mat<ldbl> v_lvmat;            // varimax回転後の因子負荷量行列
  mat<ldbl> v_avmat;            // varimax回転後の因子スコア計算用係数行列
  mat<ldbl> v_pmat;             // varimax回転用因子変換行列
  mat<ldbl> v_lpmat;            // promax回転後の因子負荷量行列
  mat<ldbl> v_mmat;             // promax回転後の因子構造行列
  mat<ldbl> v_apmat;            // promax回転後の因子スコア計算用係数行列
  mat<ldbl> v_tprmat;           // promax回転用因子変換行列
  mat<ldbl> v_ttmat;            // promax回転後の因子軸間相関行列
  int v_np;                     // パラメータ数
  int v_nq;                     // 因子数
  int v_n;                      // 例数
  int v_nrfact;                 // 因子抽出の反復計算回数
  int v_nrrot;                  // 因子軸回転の反復計算回数
  ldbl v_sllr;                  // 特殊対数尤度比
  rottype v_rtype;              // 因子軸回転の種類
  restype v_resfact;            // 因子抽出結果の種類
  restype v_resrot;             // 因子軸回転結果の種類
  void cal_fact(option& opt);   // 因子を抽出する
  ldbl init_dv(option& opt, vec<ldbl>& dv, const mat<ldbl>& rmat) const;
  // 独自性初期値を設定する
  restype fact_ml(option& opt, vec<ldbl>& dv, const mat<ldbl>& rmat);
  // 最尤法による因子の抽出
  int cal_lmat(option& opt, mat<ldbl>& qmat, int nq, ldbl lamda, const vec<ldbl>& dv, const mat<ldbl>& rmat);
  // 因子負荷行列を計算する
  restype fact_pf(option& opt, vec<ldbl>& dv, const mat<ldbl>& rmat);
  // 主因子法による因子の抽出
  ldbl fsllr(const mat<ldbl>& lmat, const mat<ldbl>& rmat, const vec<ldbl>& dv) const;
  // 特殊対数尤度比を返す
  restype varimax(int& nr, mat<ldbl>& lmat, ldbl eps, int rmax) const;
  // 因子負荷量の規準varimax回転
  void vrot(mat<ldbl>& lmat) const;
  // varimax回転主要部
  restype promax(mat<ldbl>& lmat, mat<ldbl>& tprmat, int pm) const;
  // 因子負荷量のpromax回転
  void lmat2amat(mat<ldbl>& amat) const;
  // 因子負荷量を因子スコア計算用係数(回帰係数)に変換する
public:
  varclass(size_t ni);          // コンストラクタ
  ~varclass()                   // デストラクタ
  {}
  const mat<ldbl>& smat() const // 積和行列を返す
  { return v_smat; }
  const vecldbl& mevec() const  // 平均値ベクトルを返す
  { return v_mevec; }
  const vecldbl& vrvec() const  // 分散ベクトルを返す
  { return v_vrvec; }
  const mat<ldbl>& lmat() const // 因子負荷量行列を返す
  { return v_lmat; }
  const mat<ldbl>& amat() const // 因子スコア計算用係数行列を返す
  { return v_amat; }
  const mat<ldbl>& lvmat() const // varimax回転後の因子負荷量行列を返す
  { return v_lvmat; }
  const mat<ldbl>& avmat() const // varimax回転後の因子スコア計算用係数行列を返す
  { return v_avmat; }
  const mat<ldbl>& pmat() const  // varimax回転用因子変換行列を返す
  { return v_pmat; }
  const mat<ldbl>& lpmat() const // promax回転後の因子負荷量行列を返す
  { return v_lpmat; }
  const mat<ldbl>& mmat() const  // promax回転後の因子構造行列を返す
  { return v_mmat; }
  const mat<ldbl>& apmat() const // promax回転後の因子スコア計算用係数行列を返す
  { return v_apmat; }
  const mat<ldbl>& tprmat() const // promax回転用因子変換行列を返す
  { return v_tprmat; }
  const mat<ldbl>& ttmat() const // promax回転後の因子軸間相関行列を返す
  { return v_ttmat; }
  int np() const                // パラメータ数を返す
  { return v_np; }
  int nq() const                // 因子数を返す
  { return v_nq; }
  int n() const                 // 例数を返す
  { return v_n; }
  int nrfact() const            // 因子抽出の反復計算回数を返す
  { return v_nrfact; }
  int nrrot() const             // 因子軸回転の反復計算回数を返す
  { return v_nrrot; }
  ldbl sllr() const             // 累積対数対数尤度比を返す
  { return v_sllr; }
  ldbl eval0(int idx) const     // 固有値初期値を返す
  { return (0 <= idx && idx <= nq()) ? evalvec0[idx] : 0.0; }
  ldbl eval(int idx) const      // 固有値を返す
  { return (0 <= idx && idx <= nq()) ? evalvec[idx] : 0.0; }
  ldbl l(int i, int k) const    // 因子負荷量を返す
  { return (0 <= i && i <= np() && 0 <= k && k <= nq()) ? v_lmat(i, k) : 0.0; }
  double vx(int idx) const      // 変数Xの分散を返す
  { return (0 <= idx && idx <= np()) ? v_smat(idx, idx) : 0.0; }
  rottype rtype() const         // 因子軸回転の種類を返す
  { return v_rtype; }
  restype resfact() const       // 因子抽出結果の種類を返す
  { return v_resfact; }
  restype resrot() const        // 因子軸回転結果の種類を返す
  { return v_resrot; }
  void cal(option& opt, ldata& datalst); // データを集計して各種統計量を求める
};

varclass::varclass(size_t ni) : v_smat(), v_mevec(), v_vrvec(), evalvec0(), evalvec(), v_lmat(), v_amat(), v_lvmat(), v_avmat(), v_pmat(), v_lpmat(), v_mmat(), v_apmat(), v_tprmat(), v_ttmat(), v_np(0), v_nq(0), v_n(0), v_nrfact(0), v_nrrot(0), v_sllr(0.0L), v_rtype(rotnon), v_resfact(resnon), v_resrot(resnon)
// コンストラクタ
{
  v_np = ni;
  size_t dim = v_np + 1;
  v_smat.resize(dim, dim);
  v_mevec.resize(dim, 0.0L);  v_vrvec.resize(dim, 0.0L);
  evalvec0.resize(dim, 0.0L);  evalvec.resize(dim, 0.0L);
  v_lmat.resize(dim, dim);  v_amat.resize(dim, dim);
  v_lvmat.resize(dim, dim);  v_avmat.resize(dim, dim);
  v_pmat.resize(v_np, v_np);  v_lpmat.resize(dim, dim);
  v_mmat.resize(dim, dim);  v_apmat.resize(dim, dim);
  v_tprmat.resize(v_np, v_np);  v_ttmat.resize(v_np, v_np);
}

void varclass::cal(option& opt, ldata& datalst)
// データを集計して各種統計量を求める
{
  int dim = np() + 1;
  vecdbl fdatavec(dim, 0.0);
  fdatavec[0] = 1.0;
  for (auto& cdatai : datalst)
  {
    bool blank = false;
    for (int j = 1; blank == false && j < dim; ++j)
    {
      fdatavec[j] = cdatai.fdata(j - 1);
      if (fdatavec[j] == BLANKf) blank = true;
    }
    if (blank)
    { cdatai.sdata(0, BLANKs); } // 後の計算のために最初の項目をブランクにする
    else
    {
      for (int i = 0; i < dim; ++i)
      {
        double x = fdatavec[i];
        for (int j = i; j < dim; ++j) v_smat.add(i, j, (ldbl)(x * fdatavec[j]));
      }
    }
  }
  v_n = (int)(v_smat(0, 0));
  if (v_n < 1)                  // 因子抽出不可能
  {
    v_resfact = reserr;
    return;
  }
  for (int i = 0; i < np(); ++i)
  {
    for (int j = i + 1; j < dim; ++j) v_smat.assign(j, i, v_smat(i, j));
  }
  v_smat.sweep(0);              // 単純積和行列→積和行列
  for (int i = 0; i < dim; ++i)
  {
    v_mevec[i] = v_smat(0, i);
    v_vrvec[i] = (v_n > 1) ? v_smat(i, i) / (ldbl)(v_n - 1) : 0.0L;
  }
  cal_fact(opt);
}

void varclass::cal_fact(option& opt) // 因子を抽出する
{
  mat<ldbl> rmat(v_smat.sub(1, 1, np(), np()));
  smat2rmat(rmat);
  vec<ldbl> dv(np());
  ldbl d2c = init_dv(opt, dv, rmat);
  switch (opt.meth())
  {
  case methml : v_resfact = fact_ml(opt, dv, rmat);
    break;
  default : v_resfact = fact_pf(opt, dv, rmat);
  }
  if (v_resfact == reserr) return; // 因子抽出失敗
  v_lmat.resize(np() + 1, nq() + 1);  v_amat.resize(np() + 1, nq() + 1);
  v_lvmat.resize(np() + 1, nq() + 1);  v_avmat.resize(np() + 1, nq() + 1);
  v_lpmat.resize(np() + 1, nq() + 1);  v_apmat.resize(np() + 1, nq() + 1);
  v_sllr = fsllr(v_lmat, rmat, dv);
  v_amat = v_lmat;
  lmat2amat(v_amat); // 因子負荷量行列を因子スコア計算用係数行列に変換する
  if (opt.vm() == false && opt.pm() < 1) return; // 因子軸の回転せず
  v_rtype = rotvm;  v_nrrot = 0;
  mat<ldbl> lmat0(v_lmat.sub(1, 1, np(), nq()));
  v_lvmat = v_lmat;
  v_resrot = varimax(v_nrrot, v_lvmat, EPSdef, RMAXdef);
  v_avmat = v_lvmat;
  lmat2amat(v_avmat);
  mat<ldbl> mmat(v_lvmat.sub(1, 1, np(), nq()));
  mat<ldbl> tmat(lmat0);
  mat<ldbl> llmat(tmat.trn() * lmat0);
  mat<ldbl> lmmat(tmat * mmat);
  v_pmat = llmat.inv() * lmmat;
  if (opt.pm() < 1) return;     // promax回転せず
  v_rtype = rotpm;
  v_lpmat = v_lvmat;  v_tprmat.resize(nq(), nq());
  promax(v_lpmat, v_tprmat, opt.pm());
  mat<ldbl> trmat(v_tprmat);
  v_ttmat = trmat.trn() * v_tprmat;
  v_ttmat.inv();
  mat<ldbl> fcormat(nq() + 1, nq() + 1);
  for (int i = 1; i <= nq(); ++i)
  {
    for (int j = 1; j <= nq(); ++j) fcormat.assign(i, j, v_ttmat(i - 1, j - 1));
  }
  v_mmat = v_lpmat * fcormat;
  v_apmat = v_mmat;
  lmat2amat(v_apmat);
}

ldbl varclass::init_dv(option& opt, vec<ldbl>& dv, const mat<ldbl>& rmat) const
// 独自性初期値を設定する
{
  ldbl d2c = 1.0;
  switch (opt.d2())             // 独自性初期値の種類
  {
  case d2coef : d2c = opt.d2c(); // 係数指定
    break;
  case d2q : d2c = (0 < opt.q() && opt.q() < np()) ? (ldbl)(opt.q()) : 1.0L;
    d2c = 1.0L - d2c / (ldbl)(2 * np()); // 因子数に基づく推定値
  }
  mat<ldbl> rtmat(rmat);
  rtmat.inv();
  for (int i = 0; i < np(); ++i)
  {
    ldbl di = (rtmat(i, i) > (ldbl)FDeps) ? d2c * sqrt(1.0L / rtmat(i, i)) : 0.9L;
    dv.assign(i, di);
  }
  return d2c;                   // 上限値にかける係数を返す
}

restype varclass::fact_ml(option& opt, vec<ldbl>& dv, const mat<ldbl>& rmat)
// 最尤法による因子の抽出
{
  mat<ldbl> rtmat(rmat);
  mat<ldbl> qmat(np());
  vec<ldbl> ev(np());
  ldbl lamda = 1.0L;
  if (opt.lamda() != BLANKf && opt.lamda() > 1.0) lamda = opt.lamda();
  ldbl res = 0.0L;
  v_nq = np();                  // 因子数初期値
  v_nrfact = 0;                 // 繰り返し計算数
  restype resml = resnon;
  while (resml == resnon)
  {
    ++v_nrfact;
    v_nq = cal_lmat(opt, qmat, v_nq, lamda, dv, rmat);
    if (v_nrfact == 1)
    { for (int i = 0; i <= v_nq; ++i) evalvec0[i] = evalvec[i]; }
    res = 0.0L;
    for (int i = 0; i < np(); ++i)
    {
      ldbl hi2 = 0.0L;
      for (int k = 0; k < v_nq; ++k) hi2 += v_lmat(i + 1, k + 1) * v_lmat(i + 1, k + 1);
      ldbl hi2pre = 1.0 - dv(i) * dv(i);
      ldbl er = hi2pre - hi2;
      res += er * er;
      dv.assign(i, sqrt(1.0L - hi2));
    }
    if (res < opt.eps()) resml = resok;
    else if (v_nrfact >= opt.rmax()) resml = resrmax;
  }
  return resml;                 // 因子抽出結果フラッグを返す
}

int varclass::cal_lmat(option& opt, mat<ldbl>& qmat, int nq, ldbl lamda, const vec<ldbl>& dv, const mat<ldbl>& rmat)
// 因子負荷行列を計算する
{
  mat<ldbl> rtmat(rmat);
  evalvec[0] = 0.0L;
  for (int i = 0; i < np(); ++i)
  {
    ldbl di = dv(i);
    for (int j = 0; j < np(); ++j) rtmat.assign(i, j, rmat(i, j) / di / dv(j));
    evalvec[0] += rtmat(i, i);
  }
  vec<ldbl> ev(np());
  int k = 0;
  if (0 < opt.q() && opt.q() <= np()) nq = opt.q(); // 最大因子数指定の時
  while (k < nq && (evalvec[k + 1] = eigenud(rtmat, ev)) > lamda)
  {
    for (int j = 0; j < np(); ++j) qmat.assign(j, k, ev(j));
    ++k;
  }
  if (k < 1) return 0;      // 因子抽出できず
  nq = k;                   // level指定の時、因子数は初期値以下になる
  vecldbl lamvec(nq, 0.0);
  for (k = 0; k < nq; ++k) lamvec[k] = sqrt(evalvec[k + 1] - 1.0);
  for (int i = 0; i < np(); ++i)
  {
    ldbl di = dv(i);
    for (k = 0; k < nq; ++k) v_lmat.assign(i + 1, k + 1, di * qmat(i, k) * lamvec[k]);
  }
  return nq;                    // 抽出した固有値の数を返す
}

restype varclass::fact_pf(option& opt, vec<ldbl>& dv, const mat<ldbl>& rmat)
// 主因子法による因子の抽出
{
  v_nq = np();                  // 因子数初期値
  ldbl lamda = 1.0L;
  vec<ldbl> ev(np());
  if (0 < opt.q() && opt.q() <= np()) { v_nq = opt.q();  lamda = 0.0L; }
  else if (opt.lamda() != BLANKf) lamda = opt.lamda();
  evalvec0[0] = evalvec[0] = (ldbl)np();
  mat<ldbl> rtmat(rmat);
  int k = 0;
  while (k < v_nq && (evalvec0[k + 1] = eigenud(rtmat, ev)) > lamda) ++k;
  v_nq = k;                     // 因子数決定
  if (v_nq < 1) return reserr;  // 因子抽出できず
  mat<ldbl> rlmat(rmat);
  if (opt.rmax() > 0)
  {
    for (int i = 0; i < np(); ++i) rlmat.assign(i, i, 1.0 - dv(i) * dv(i));
  }
  ldbl res = 0.0L;
  restype respf = resnon;
  v_nrfact = 0;                 // 繰り返し計算数
  while (respf == resnon)
  {
    ++v_nrfact;
    rtmat = rlmat;
    for (int k = 1; k <= v_nq; ++k)
    {
      evalvec[k] = eigenud(rtmat, ev);
      ldbl w = sqrt(evalvec[k]);
      for (int i = 1; i <= np(); ++i) v_lmat.assign(i, k, ev(i - 1) * w);
    }
    res = 0.0L;
    for (int i = 1; i <= np(); ++i)
    {
      ldbl hi2 = 0.0L;
      for (int k = 1; k <= v_nq; ++k) hi2 += v_lmat(i, k) * v_lmat(i, k);
      dv.assign(i - 1, (hi2 < 1.0) ? sqrt(1.0 - hi2) : 0.0L);
      ldbl er = rlmat(i - 1, i - 1) - hi2;
      res += er * er;
      rlmat.assign(i - 1, i - 1, hi2);
    }
    if (res < opt.eps()) respf = resok;
    else if (v_nrfact >= opt.rmax()) respf = resrmax;
  }
  return respf;                 // 因子抽出結果フラッグを返す
}

ldbl varclass::fsllr(const mat<ldbl>& lmat, const mat<ldbl>& rmat, const vec<ldbl>& dv) const
// 特殊対数尤度比を返す
{
  int np = lmat.r() - 1, nq = lmat.c() - 1;
  mat<ldbl> ltmat(lmat.sub(1, 1, np, nq));
  mat<ldbl> lttmat(ltmat);
  lttmat.trn();
  mat<ldbl> rlmat(ltmat * lttmat);
  for (int i = 0; i < np; ++i) rlmat.assign(i, i, rlmat(i, i) + dv(i) * dv(i));
  ldbl det = rmat.det();
  ldbl lnrr = (fabs(det) > (ldbl)FDeps) ? log(rlmat.det() / det) : 0.0L;
  rlmat.inv();
  mat<ldbl> tmat(rlmat * rmat);
  ldbl tl = 0.0;
  for (int i = 0; i < np; ++i) tl += tmat(i, i);
  return lnrr + tl - (ldbl)np; // 特殊対数尤度を返す
}

restype varclass::varimax(int& nr, mat<ldbl>& lmat, ldbl eps, int rmax) const
// 因子負荷量の規準varimax回転
{
  nr = 0;
  int np = (int)(lmat.r()) - 1, nq = (int)(lmat.c()) - 1; // 変数の数と因子数
  if (np < 2 || nq < 2) return reserr;   // 回転不可能
  mat<ldbl> amat(lmat.sub(1, 1, np, nq));
  for (int i = 0; i < np; ++i)
  {
    ldbl norm = 0.0;
    for (int j = 0; j < nq; ++j) norm += amat(i, j) * amat(i, j);
    norm = sqrt(norm);
    lmat.assign(i + 1, 0, norm);
    for (int j = 0; j < nq; ++j) amat.assign(i, j, (norm > 0.0) ? amat(i, j) / norm : 0.0);
  }
  ldbl sign = 0.0L, sigo = 0.0L, sigm = 0.0L, sigm2 = 0.0L, m = 0.0L;
  restype resvm = resnon;
  while (resvm == resnon)
  {
    ++nr;
    sigo = sign;
    vrot(amat);
    sign = 0.0L;
    for (int j = 0; j < nq; ++j)
    {
      sigm = sigm2 = 0.0L;
      for (int i = 0; i < np; ++i)
      {
        m = amat(i, j) * amat(i, j);
        sigm += m;  sigm2 += m * m;
      }
      sign += sigm2 - sigm * sigm / (ldbl)np;
      if (fabs(sign - sigo) < eps * (fabs(sign) + fabs(sigo))) resvm = resok;
      else if (nr >= rmax) resvm = resrmax;
    }
  }
  for (int i = 1; i <= np; ++i)
  {
    ldbl norm = lmat(i, 0);
    for (int j = 1; j <= nq; ++j) lmat.assign(i, j, amat(i - 1, j - 1) * norm);
  }
  return resvm;                 // 回転成功
}

void varclass::vrot(mat<ldbl>& lmat) const
// varimax回転主要部
{
  int np = (int)(lmat.r()), nq = (int)(lmat.c()); // 変数の数と因子数
  for (int j = 0; j < nq - 1; ++j)
  {
    for (int k = j + 1; k < nq; ++k)
    {
      ldbl a = 0.0, b = 0.0, c = 0.0, d = 0.0;
      for (int i = 0; i < np; ++i)
      {
        ldbl lj = lmat(i, j), lk = lmat(i, k);
        ldbl ljjkk = lj * lj - lk * lk, ljk = lj * lk;
        a += ljjkk;  b += ljk;
        c += ljjkk * ljjkk - 4.0 * ljk * ljk;  d += ljjkk * ljk;
      }
      b *= 2.0;  d *= 4.0;
      d -= 2.0 * a * b / (ldbl)np;
      b = c - (a * a - b * b) / (ldbl)np;
      ldbl h = atan2(d, b) / 4.0;
      for (int i = 0; i < np; ++i)
      {
        ldbl lj = lmat(i, j), lk = lmat(i, k);
        lmat.assign(i, j, lj * cos(h) + lk * sin(h));
        lmat.assign(i, k, lk * cos(h) - lj * sin(h));
      }
    }
  }
}

restype varclass::promax(mat<ldbl>& lmat, mat<ldbl>& tprmat, int pm) const
// 因子負荷量のpromax回転
{
  int np = (int)(lmat.r()) - 1, nq = (int)(lmat.c()) - 1; // 変数の数と因子数
  if (np < 2 || nq < 2) return reserr;
  mat<ldbl> amat(lmat.sub(1, 1, np, nq));
  mat<ldbl> cmat(np, nq);       // ターゲット行列
  ldbl pmf = (ldbl)pm;
  for (int j = 0; j < np; ++j)
  {
    ldbl norm = 0.0L;
    for (int k = 0; k < nq; ++k) norm += amat(j, k) * amat(j, k);
    norm = sqrt(norm);
    for (int k = 0; k < nq; ++k)
    {
      ldbl ljk = (norm > 0.0) ? amat(j, k) / norm : 0.0L;
      cmat.assign(j, k, sgn(ljk) * pow(fabs(ljk), pmf));
    }
  }
  mat<ldbl> tmat(amat);
  mat<ldbl> aamat(tmat.trn() * amat);
  aamat.inv();
  mat<ldbl> acmat(tmat * cmat);
  mat<ldbl> trmat(aamat * acmat);
  mat<ldbl> ttmat(trmat);
  mat<ldbl> dmat(ttmat.trn() * trmat);
  dmat.inv();
  for (size_t j = 0; j < dmat.c(); ++j)
  {
    for (size_t i = 0; i < dmat.r(); ++i)
    {
      if (i == j) dmat.assign(i, j, (dmat(i, j) > 0.0) ? sqrt(dmat(i, j)) : 0.0L);
      else dmat.assign(i, j, 0.0L);
    }
  }
  tprmat = trmat * dmat;
  tmat.trn();
  tmat = amat * tprmat;
  for (int j = 1; j <= np; ++j)
  { for (int k = 1; k <= nq; ++k) lmat.assign(j, k, tmat(j - 1, k - 1)); }
  return resok;                 // 回転成功
}

void varclass::lmat2amat(mat<ldbl>& amat) const
// 因子負荷量を因子スコア計算用係数(回帰係数)に変換する
{
  int np = (int)(amat.r()) - 1, nq = (int)(amat.c()) - 1;
  mat<ldbl> rmat(v_smat.sub(1, 1, np, np));
  smat2rmat(rmat);
  rmat.inv();
  mat<ldbl> tmat(rmat * amat.sub(1, 1, np, nq));
  vecldbl wvec(np, 0.0);
  int df = n() - 1;
  for (int i = 0; i < np; ++i)
  {
    ldbl ss = v_smat(i + 1, i + 1);
    wvec[i] = (ss > 0.0L) ? sqrt((ldbl)df / ss) : 0.0L;
  }
  for (int j = 1; j <= nq; ++j)
  {
    amat.assign(0, j, 0.0L);
    for (int i = 1; i <= np; ++i)
    {
      amat.assign(i, j, tmat(i - 1, j - 1) * wvec[i - 1]);
      amat.assign(0, j, amat(0, j) - amat(i, j) * v_smat(0, i));
    }
  }
}

class factclass                 // factメインクラス
{
private:
  string fname;                 // ファイル名
  idd id;                       // データ識別情報
  ldata datalst;                // データリスト
  void put_resid(ostream& fout);
  // 解析結果のタイトルとデータ識別情報を出力する
  void put_factor(ostream& fout, option& opt, const varclass& var);
  // 因子分析結果を出力する
  void put_eval(ostream& fout, const varclass& var) const;
  // 固有値を出力する。
  void put_fload(ostream& fout, option& opt, const string& mtitle, const mat<ldbl>& lmat, bool fl) const;
  // 因子負荷量行列を表形式で出力する。
  void put_rfact(ostream& fout, const string& mtitle, const mat<ldbl>& amat0, const varclass& var) const;
  // 因子スコア間相関係数とCronbachのα係数を出力する。
  void put_rot(ostream& fout, option& opt, const varclass& var) const;
  // 回転後の因子分析結果を出力する
  void put_lt(ostream& fout, const varclass& var) const;
  // 尤度比検定結果を出力する。
  void sort_lvecvec(vector<vecldbl >& lvecvec, vector<vecldbl*>& lvecpvec) const;
  // 因子負荷量ベクトルを因子順に並べ替える(Schellソート)
  int cmp_vec(const vecldbl* lvecp1, const vecldbl* lvecp2) const;
  // 因子負荷量ベクトルの比較
  void put_pmat(ostream& fout, const string& mtitle, const mat<ldbl>& pmat) const;
  // 回転用因子変換行列を出力する。
  void put_dns(option& opt, const varclass& var);
  // 因子スコアをデータファイル形式で出力する
  void plot_flplot(option& opt, const varclass& var);
  // Ngraph用因子負荷量のプロットデータを出力する
public:
  factclass() : fname(), id(), datalst() // コンストラクタ
  {}
  ~factclass()                  // デストラクタ
  {}
  int get(const char* c_fname, option& opt); // データを読み込み保存する
  void put(option& opt);                     // 集計結果を出力する
};

void factclass::put_resid(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) << "\n";
  }
}

void factclass::put_factor(ostream& fout, option& opt, const varclass& var)
// 因子分析結果を出力する
{
  if (opt.res("fa"))
  {
    fout << "\n・因子分析結果:";
    switch (opt.meth())
    {
    case methml : fout << "最尤法  反復回数=" << var.nrfact();
      break;
    default :
      if (opt.rmax() < 1) fout << "主成分法";
      else fout << "主因子法  反復回数=" << var.nrfact();
    }
    if (opt.rmax() > 0 && var.resfact() == resrmax) fout << ":最大回数";
    fout << "\n  因子抽出条件:";
    if (0 < opt.q() && opt.q() <= var.np()) fout << "最大因子数=" << opt.q();
    else fout << "初期固有値＞" << opt.lamda();
    fout << "  独自性初期値=";
    switch (opt.d2())
    {
    case d2coef : fout << "上限値×" << opt.d2c();
      break;
    case d2q : fout << "因子数に基づく推定値";
      break;
    default : fout << "上限値";
    }
    fout << "\n";
    put_eval(fout, var);
    if (opt.meth() == methml && opt.res("lt")) put_lt(fout, var);
  }
  if (opt.res("fl")) put_fload(fout, opt, "・因子負荷量", var.lmat(), true);
  if (opt.res("fc")) putcoef(fout, "・因子スコア計算用係数(回帰係数)", "因子", "x", var.amat());
  if (opt.res("fr")) put_rfact(fout, "・因子スコア間相関係数とCronbachのα係数", var.amat(), var);
}
  
void factclass::put_eval(ostream& fout, const varclass& var) const
// 固有値を出力する。
{
  fout << "因子  初期固有値  収束固有値    負荷量SS  累積負荷量      寄与率  累積寄与率\n";
  constexpr int LENline = 76;   // 区切り線の長さ
  putline(fout, LENline);
  ldbl slss = 0.0L;
  for (int k = 1; k <= var.nq(); ++k)
  {
    fout << "F" << setw(3) << k << setw(12) << var.eval0(k);
    fout << setw(12) << var.eval(k);
    ldbl lss = 0.0L;
    for (int i = 1; i <= var.np(); ++i) lss += var.l(i, k) * var.l(i, k);
    slss += lss;
    fout << setw(12) << lss << setw(12) << slss;
    fout << setw(12) << lss / (ldbl)(var.np());
    fout << setw(12) << slss / (ldbl)(var.np()) << "\n";
  }
  putline(fout, LENline);
}

void factclass::put_lt(ostream& fout, const varclass& var) const
// 尤度比検定結果を出力する。
{
  ldbl llr = -0.5L * (ldbl)(var.n()) * var.sllr();
  fout << "・因子数に関するズレの検定  尤度比:LR=" << exp(llr) << "  ln(LR)=" << llr << "\n";
  ldbl csqr = ((ldbl)var.n() - (ldbl)(2 * var.np() + 5) / 6.0L - 2.0L * (ldbl)(var.nq()) / 3.0L) * var.sllr();
  int df = ((var.np() - var.nq()) * (var.np() - var.nq()) - (var.np() + var.nq())) / 2;
  if (sgn(csqr) < 0 || df < 1)
  { fout << "  検定不可能！\n"; }
  else
  {
    ldbl pr = ctop(csqr, df);
    fout << "  χ2乗=" << csqr << "  自由度=" << df << "  有意確率p=" << pr << pind(pr) << "\n";
  }
}

void factclass::put_fload(ostream& fout, option& opt, const string& mtitle, const mat<ldbl>& lmat, bool fl) const
// 因子負荷量行列を表形式で出力する。
{
  int np = lmat.r() - 1, nq = lmat.c() - 1;
  vecldbl lvec(nq + 3, 0.0);
  vector<vecldbl> lvecvec(np, lvec); // ソート用因子負荷量ベクトル
  vector<vecldbl*> lvecpvec(np); // ソート用ポインターベクトル
  for (int i = 0; i < np; ++i)   // 因子負荷量に変数番号と因子番号を付加する
  {
    ldbl lmax = 0.0L;
    int jmax = 0;
    (lvecvec[i])[2] = 0.0L;
    for (int j = 3; j < nq + 3; ++j)
    {
      (lvecvec[i])[j] = lmat(i + 1, j - 2);
      ldbl labs = fabs(lmat(i + 1, j - 2));
      if (labs > lmax) { lmax = labs;  jmax = j; }
    }
    (lvecvec[i])[0] = (ldbl)i;
    (lvecvec[i])[1] = (ldbl)jmax;
    lvecpvec[i] = &(lvecvec[i]);
  }
  if (opt.sort()) sort_lvecvec(lvecvec, lvecpvec);
  vecldbl l2vec(nq + 1, 0.0L);
  for (int j = 0; j <= nq; ++j) l2vec[j] = 0.0L;
  constexpr int Cmax = 6;       // 行列表横の最大数
  int low = 1, up = 0, eopl = (fl) ? 2 : 0;
  bool eop = false;             // 表の最後出力フラッグ
  fout << "\n" << mtitle;
  int prec0 = fout.precision(); // 現在の浮動小数点精度を保存する
  while (eop == false)
  {
    int nc = low + Cmax - 1;
    if (low > nq) eop = true;
    else if (nc >= nq + eopl) { up = nq;  eop = true; }
    else if (nc >= nq) up = nq;
    else up = nc;
    fout << "\n変数";
    for (int j = low; j <= up; j++) fout << "      因子" << setw(2) << j;
    if (eop && fl) fout << "      共通性      独自性";
    fout << "\n";
    int llen = 4 + 12 * (up - low + 1 + ((eop) ? eopl : 0));
    putline(fout, llen);
    for (int i = 0; i < np; ++i)
    {
      vecldbl* lvecpi = lvecpvec[i];
      fout << "x" << setw(3) << (*lvecpi)[0] + 1;
      fout.precision(6);  fout.setf(ios::fixed);
      for (int j = low; j <= up; ++j)
      {
        fout << setw(12) << (*lvecpi)[j + 2];
        ldbl l2 = (*lvecpi)[j + 2] * (*lvecpi)[j + 2];
        (*lvecpi)[2] += l2;
        l2vec[j] += l2;
      }
      if (eop && fl) fout << setw(12) << (*lvecpi)[2] << setw(12) << 1.0L - (*lvecpi)[2];
      fout.precision(prec0);  fout.unsetf(ios::fixed);
      fout << "\n";
    }
    putline(fout, llen);
    if (fl)
    {
      fout << "寄与率";
      fout.precision(6);  fout.setf(ios::fixed);
      ldbl r2 = 0.0L;
      for (int j = low; j <= up; ++j)
      {
        r2 = l2vec[j] / (ldbl)np;
        if (j == low) fout << setw(10) << r2;
        else fout << setw(12) << r2;
        l2vec[0] += l2vec[j];
      }
      if (eop)
      {
        r2 = l2vec[0] / (ldbl)np;
        if (low > nq) fout << setw(10) << r2;
        else fout << setw(12) << r2;
        fout << setw(12) << 1.0L - r2;
      }
      fout.precision(prec0);  fout.unsetf(ios::fixed);
      fout << "\n";
    }
    low = up + 1;
  }
}

void factclass::put_rfact(ostream& fout, const string& mtitle, const mat<ldbl>& amat0, const varclass& var) const
// 因子スコア間相関係数とCronbachのα係数を出力する。
{
  constexpr int Cmax = 6;       // 行列表横の最大個数
  mat<ldbl> amat(amat0);
  int np = amat.r() - 1, nq = amat.c() - 1, low = 1, up = 0, eopl = 1;
  ldbl df = (ldbl)(var.n() - 1);
  for (int i = 0; i <= nq; ++i)
  { amat.assign(i, 0, 0.0L);  amat.assign(0, i, 0.0L); }
  for (int i = nq + 1; i <= np; ++i) amat.assign(i, 0, 0.0L);
  mat<ldbl> atmat(amat);
  mat<ldbl> attmat(atmat.trn() * var.smat());
  mat<ldbl> rmat(attmat * amat);
  bool eop = false;             // 表の最後の出力フラッグ
  fout << "\n" << mtitle;
  int prec0 = fout.precision(); // 現在の浮動小数点精度を保存する
  while (eop == false)
  {
    int nc = low + Cmax - 1;
    if (low > nq) eop = true;
    else if (nc >= nq + eopl) { up = nq;  eop = true; }
    else if (nc >= nq) up = nq;
    else up = nc;
    fout << "\n      ";
    for (int j = low; j <= up; j++) fout << "      因子" << setw(2) << j;
    if (eop) fout << "        分散";
    fout <<  "\n";
    int llen = 6 + 12 * (up - low + 1 + ((eop) ? eopl : 0));
    putline(fout, llen);
    for (int i = 1; i <= nq; ++i)
    {
      fout << "因子" << setw(2) << i;
      fout.precision(6);  fout.setf(ios::fixed);
      for (int j = low; j <= up; ++j)
      {
        if (i == j)
        { fout << setw(12) << 1.0; }
        else
        {
          ldbl r = rmat(i, i) * rmat(j, j);
          r = (r > 0.0L) ? rmat(i, j) / sqrt(r) : 0.0L;
          fout << setw(12) << r;
        }
      }
      if (eop) fout << setw(12) << rmat(i, i) / df;
      fout.precision(prec0);  fout.unsetf(ios::fixed);
      fout << "\n";
    }
    putline(fout, llen);
    fout << "α係数";
    fout.precision(6);  fout.setf(ios::fixed);
    for (int j = low; j <= up; ++j)
    {
      ldbl sxx = 0.0;
      for (int i = 1; i <= np; ++i) sxx += amat(i, j) * var.vx(i) * amat(i, j);
      ldbl alpha = (rmat(j, j) > 0.0L) ? (ldbl)np / (ldbl)(np - 1) * (1.0L - sxx / rmat(j, j)) : 0.0L;
      fout << setw(12) << alpha;
    }
    fout.precision(prec0);  fout.unsetf(ios::fixed);
    fout << "\n";
    low = up + 1;
  }
}

void factclass::put_rot(ostream& fout, option& opt, const varclass& var) const
// 回転後の因子分析結果を出力する
{
  ostringstream ostr;
  if (opt.res("fl"))
  {
    ostr << "・規準varimax回転後の因子負荷量  反復回数=" << var.nrrot();
    if (var.resrot() == reserr) ostr << ":回転不可能！";
    else if (var.resrot() == resrmax) ostr << ":最大回数";
    put_fload(fout, opt, ostr.str(), var.lvmat(), true);
  }
  if (var.resrot() == reserr) return;
  if (opt.res("fc")) putcoef(fout, "・規準varimax回転後の因子スコア計算用係数(回帰係数)", "因子", "x", var.avmat());
  if (opt.res("ft")) put_pmat(fout, "・規準varimax回転用因子変換行列", var.pmat());
  if (opt.res("fr")) put_rfact(fout, "・規準varimax回転後の因子スコア間相関係数とCronbachのα係数", var.avmat(), var);
  if (opt.pm() < 1) return; // promax回転せず
  ostr.str("");
  ostr << "・promax回転後の因子負荷量(ベキ乗数=" << opt.pm() << ")";
  if (opt.res("fl")) put_fload(fout, opt, ostr.str(), var.lpmat(), false);
  if (opt.res("fc"))
  {
    put_fload(fout, opt, "・promax回転後の因子構造行列(回帰係数計算用)", var.mmat(), false);
    putcoef(fout, "・promax回転後の因子スコア計算用係数(回帰係数)", "因子", "x", var.apmat());
  }
  if (opt.res("ft")) put_pmat(fout, "・promax回転用因子変換行列", var.tprmat());
  if (opt.res("fr"))
  {
    put_pmat(fout, "・promax回転後の因子軸間相関行列(r=cos(θ))", var.ttmat());
    put_rfact(fout, "・promax回転後の因子スコア間相関係数とCronbachのα係数", var.apmat(), var);
  }
}

void factclass::sort_lvecvec(vector<vecldbl >& lvecvec, vector<vecldbl*>& lvecpvec) const
// 因子負荷量ベクトルを因子順に並べ替える(Schellソート)
{
  int np = lvecvec.size();
  for (int gap = np / 2; gap > 0; gap /=2)
  {
    for (int i = gap; i < np; ++i)
    {
      vecldbl* lvecpi = lvecpvec[i];
      int j = 0;
      for (j = i - gap; j >= 0 && cmp_vec(lvecpi, lvecpvec[j]) < 0; j -= gap) lvecpvec[j + gap] = lvecpvec[j];
      lvecpvec[j + gap] = lvecpi;
    }
  }
}

int factclass::cmp_vec(const vecldbl* lvecp1, const vecldbl* lvecp2) const
// 因子負荷量ベクトルの比較
{
  int jmax1 = (int)((*lvecp1)[1]), jmax2 = (int)((*lvecp2)[1]);
  if (jmax1 < jmax2) return -1; // 因子番号は小さい順
  if (jmax1 > jmax2) return 1;
  ldbl labs1 = fabs((*lvecp1)[jmax1]), labs2 = fabs((*lvecp2)[jmax2]);
  if (labs1 > labs2) return -1; // 因子負荷量は絶対値の大きい順
  if (labs1 < labs2) return 1;
  return 0;
}

void factclass::put_pmat(ostream& fout, const string& mtitle, const mat<ldbl>& pmat) const
// 回転用因子変換行列を出力する。
{
  constexpr int Cmax = 6;       // 行列表横の最大個数
  int p = pmat.r(), q = pmat.c(), low = 0, up = 0;
  fout << "\n" << mtitle;
  int prec0 = fout.precision(); // 現在の浮動小数点精度を保存する
  while (low < q)
  {
    up = (low + Cmax < q) ? low + Cmax : q;
    fout << "\n      ";
    for (int j = low; j < up; j++) fout << "      因子" << setw(2) << j + 1;
    fout <<  "\n";
    int llen = 6 + 12 * (up - low);
    putline(fout, llen);
    for (int i = 0; i < q; ++i)
    {
      fout << "因子" << setw(2) << i + 1;
      fout.precision(6);  fout.setf(ios::fixed);
      for (int j = low; j < up; ++j) fout << setw(12) << pmat(i, j);
      fout.precision(prec0);  fout.unsetf(ios::fixed);
      fout << "\n";
    }
    putline(fout, llen);
    low = up;
  }
}

void factclass::put_dns(option& opt, const varclass& var)
// 因子スコアをデータファイル形式で出力する
{
  ofstream fout;
  if (ofopen(fout, (opt.dnsname()).c_str()) == false) return; // ファイルオープン失敗
  string tstr("");
  mat<ldbl> amat(var.amat());
  if (var.rtype() == rotvm)
  {
    tstr = "(varimax回転後)";
    amat = var.avmat();
  }
  else if (var.rtype() == rotpm)
  {
    tstr = "(promax回転後)";
    amat = var.apmat();
  }
  id.put_title(fout, formtag);
  id.put_stra(fout, formtag);
  idd tid(id.title(), var.nq());
  ostringstream ostr;
  for (int i = 0; i < var.nq(); ++i)
  {
    ostr << "因子スコアF" << i + 1 << tstr;
    tid.assign_item(i, ostr.str());
    ostr.str("");
  }
  tid.put_item(fout, formtag);
  datac score(var.nq());
  for (auto& cdatai : datalst)
  {
    if (cdatai.sdata(0) == BLANKs) // 欠損あり
    { score = BLANKs; }
    else
    {
      for (int i = 1; i <= var.nq(); ++i)
      {
        ldbl z = amat(0, i);
        for (size_t j = 0; j < id.ni(); ++j) z += amat(j + 1, i) * (ldbl)(cdatai.fdata(j));
        score.assign(i - 1, z);
      }
    }
    score.put(fout);
  }
}

void factclass::plot_flplot(option& opt, const varclass& var)
// Ngraph用因子負荷量のプロットデータを出力する
{
  ofstream fout;
  if (ofopen(fout, (opt.ngpname()).c_str()) == false) return; // ファイルオープン失敗
  fout << "-s1 -zx-1,1,0.5 -zy-1,1,0.5\n";
  fout << "%";  puttitle(fout, Title);
  fout << "%\n";
  id.put_title(fout, formngp);
  id.put_stra(fout, formngp);
  for (size_t i = 0; i < id.ni(); ++i)
  {
    fout << "%変数x" << setw(3) << i + 1;
    fout << ":" << id.item(i) << "\n";
  }
  fout << "%円を描く時:-dline、プロットを描く時:-dmark,00 -o300\n";
  constexpr uint RDdef = 60;    // 座標軸の分割数
  ldbl dh = 2.0L * M_PI / (ldbl)RDdef;
  ldbl h = 0.0L;
  for (uint i = 0; i <= RDdef; ++i, h += dh) fout << cos(h) << " " <<  sin(h) << "\n";
  string tstr("");
  mat<ldbl> lmat(var.lmat());
  if (var.rtype() == rotvm)
  {
    tstr = " (varimax回転後)";
    lmat = var.lvmat();
  }
  else if (var.rtype() == rotpm)
  {
    tstr = " (promax回転後)";
    lmat = var.lpmat();
  }
  for (int j = 1; j < var.nq(); ++j)
  {
    for (int k = j + 1; k <= var.nq(); ++k)
    {
      fout << NGPdelm << "\n";
      fout << "%x:F" << j << "  y:F" << k << tstr << "\n";
      for (int i = 1; i <= var.np(); ++i) fout << lmat(i, j) << " " << lmat(i, k) << "\n";
    }
  }
}

int factclass::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 factclass::put(option& opt) // 集計結果を出力する
{
  varclass var(id.ni());        // 変数オブジェクト
  var.cal(opt, datalst);        // 各種の計算実行
  if (opt.res("res"))
  {
    put_resid(cout);
    if (opt.res("stat"))
    {
      putmsdvec(cout, var.n(), var.mevec(), var.vrvec());
      putmat(cout, "・相関行列", matr, var.n(), var.smat());
    }
    if (var.resfact() == reserr)
    {
      cout << "※因子抽出不可能！\n\n" << endl;
      return;
    }
    put_factor(cout, opt, var);
    if (var.nq() > 1 && (opt.vm() || opt.pm() > 0)) put_rot(cout, opt, var);
    cout << "\n" << endl;
  }
  if (opt.dnsname() != "") put_dns(opt, var);
  if (opt.ngpname() != "") plot_flplot(opt, var);
}

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