// tcoreg.C Version 7.44  Copyleft (c) 杉本典夫
// Time-stamp: <2025-10-29 14:59:56 tac>
//
// 時期別相関係数と回帰直線、連動係数と連動回帰直線を計算する。
//
// 使用ライブラリ:libdans.so, libfun.so
// コンパイル法  :g++ -O -s -W tcoreg tcoreg.C
//                -I/usr/local/include -L/usr/local/lib -ldans -lfun

#include <dans.H>               // dans用ヘッダファイル(fun.H含む)

using namespace dans;           // dans名前空間使用(名前空間matvec含む)

const vecstr USAGE =            // プログラム使用法表示用
{
  "tcoreg  Version 7.44         Copyleft(Vo\\oV) 2025 杉本典夫",
  "書式:tcoreg  [<オプション> …]  [<データファイル> …]",
  "機能:時期別相関係数・回帰直線と連動係数・連動回帰直線を計算する。",
  "<オプション>",
  "  -cc[=<信頼係数>]:信頼係数を指定する。(0≦信頼係数<1)",
  "    -cc:-cc=0.95(オプション省略時)",
  "  -dxrange=<最小値>-<最小値>:X変化量の最小値と最大値を指定する。",
  "  -f=<ファイル名>:指定されたファイルからオプションを読み込む。",
  "  -h[=<コマンド名>]:オンラインマニュアルを表示します。(-h:-h=man)",
  "  -ngp={ce|data|ddata|line|lline}[,…]:Ngraph用出力データを指定する。",
  "    ce:信頼・許容・予測楕円  data:実測値生データ  ddata:変化量生データ",
  "    line:回帰直線  lline:連動回帰直線  データオプション省略時:全て指定",
  "  -ngpname[=<出力ファイル名>]:Ngraph用データを出力するファイル名を指定する。",
  "    オプション省略時:Ngraph用データ出力せず  -ngpname:-ngpname=tcoreg_ngp.dat",
  "  -res=[-][{def|all|l|r|stat}][,…]:解析結果の出力項目を指定する。",
  "    -:-を付けた項目を出力しない(-res=-:全ての項目を出力しない)",
  "    def:オプション省略時と同じ項目を出力する  all:全ての項目を出力する",
  "    l:連動係数と連動回帰直線  r:相関係数と回帰直線  stat:基礎統計量",
  "    オプション省略時:自動判定",
  "  -ti[={-|+}]:時期別解析結果の出力指定。(-ti:ti=+)",
  "  -tt[={-|+}]:時期合計解析結果の出力指定。(-tt:tt=+)",
  "  -x=[<開始項目番号>][-[<終了項目番号>]]:説明変数の項目番号を指定する。",
  "  -xrange=<最小値>-<最小値>:X実測値の最小値と最大値を指定する。",
  "  -y=<開始項目番号>:目的変数の開始項目番号を指定する。"
};

const string Title = "=== 時期別相関と連動係数 ==="; // 表題
const string NGPNAMEdef = "tcoreg_ngp.dat"; // 省略時Ngraphファイル名

constexpr int LENline = 75;     // 区切り線の長さ

class option                    // オプション集合体
{
  string v_path;                // プログラムのパス
  string v_pname;               // プログラム名
  double v_cc;                  // 信頼係数
  pairdbl dxrange;              // x変化量の最小値と最大値
  string v_ngpname;             // Ngraphファイル名
  lstr ngplst;                  // Ngraph出力リスト
  mapsb ngpmap;                 // Ngraph出力の有無マップ
  lstr reslst;                  // 分析結果出力リスト
  mapsb resmap;                 // 分析結果出力の有無マップ
  bool v_ti;                    // 時期別集計結果出力フラッグ
  bool v_tt;                    // 時期合計集計結果出力フラッグ
  int v_nt;                     // 時期数
  pairdbl xrange;               // x実測値の最小値と最大値
  int v_xs;                     // 変数Xの開始項目番号
  int v_ys;                     // 変数Yの開始項目番号
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; }
  double dxmin() const          // x変化量の最小値を返す
  { return dxrange.first; }
  double dxmax() const          // x変化量の最大値を返す
  { return dxrange.second; }
  const string& ngpname() const // Ngraphデータファイル名を返す
  { return v_ngpname; }
  bool ngp(const string& code); // Ngraph出力フラッグを返す
  bool res(const string& code); // 分析結果出力フラッグを返す
  bool ti() const               // 時期別集計結果出力フラッグを返す
  { return v_ti; }
  bool tt() const               // 時期合計集計結果出力フラッグを返す
  { return v_tt; }
  int nt() const                // 時期数を返す
  { return v_nt; }
  double xmin() const           // x実測値の最小値を返す
  { return xrange.first; }
  double xmax() const           // x実測値の最大値を返す
  { return xrange.second; }
  int xs() const                // 変数Xの開始項目番号を返す
  { return v_xs; }
  int ys() const                // 変数Yの開始項目番号を返す
  { return v_ys; }
  void setdef()                 // 省略時オプション指定を設定する
  { getoptdef(*this); }
  void set(const string& optstr, char devc = OPTdev);
  // オプション指定によってオプションを設定する
  int check(size_t ni);         // オプションのチェックとセット
  int set_ngp();                // -ngpオプションのセット
  int set_res();                // -resオプションのセット
};

option::option(const char* c_pname) : v_path(), v_pname(), v_cc(CCdef), dxrange(BLANKf, BLANKf), v_ngpname(), ngplst(), ngpmap(), reslst(), resmap(), v_ti(true), v_tt(true), v_nt(1), xrange(BLANKf, BLANKf), v_xs(1), v_ys(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::ngp(const string& code) // Ngraph出力フラッグを返す
{
  mapsb_citr ngppos = ngpmap.find(code);
  if (ngppos == ngpmap.end()) return false; // code無し！
  return ngpmap[code];                      // true/falseを返す
}

bool option::res(const string& code) // 分析結果出力フラッグを返す
{
  mapsb_citr respos = resmap.find(code);
  if (respos == 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] == "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] == "dxrange") // x変化量の最小値と最大値の指定
  {
    if (getrange(dxrange.first, dxrange.second, opt[1]) == 0) // 範囲指定無
    { err = badopt(pname(), optstr, USAGE, -FDmax, FDmax); }
    else if (dxrange.first >= dxrange.second) // 範囲指定不適切
    { err = badopt(pname(), optstr, USAGE, -FDmax, FDmax); }
  }
  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] == "ngp")     // Ngraph出力項目指定
  {
    if (strsplit(ngplst, opt[1], CHRdelm) < 1) 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] == "res")     // 分析結果の出力指定
  {
    if (strsplit(reslst, opt[1], CHRdelm) < 1) err = badopt(pname(), optstr, USAGE);
  }
  else if (opt[0] == "ti")      // 時期別集計結果の出力指定
  {
    if (opt[1].empty() || opt[1] == "+") v_ti = true;
    else if (opt[1] == "-") v_ti = false;
    else err = badopt(pname(), optstr, USAGE);
  }
  else if (opt[0] == "tt")      // 時期合計集計結果の出力指定
  {
    if (opt[1].empty() || opt[1] == "+") v_tt = true;
    else if (opt[1] == "-") v_tt = false;
    else err = badopt(pname(), optstr, USAGE);
  }
  else if (opt[0] == "x")       // 説明変数の項目指定
  {
    int xe = 0;
    if (getrange(v_xs, xe, opt[1]) < 1) err = badopt(pname(), optstr, USAGE);
    else if (v_xs < 1 || v_xs > xe) err = badopt(pname(), optstr, USAGE);
    else v_nt = xe - v_xs + 1;
  }
  else if (opt[0] == "xrange")  // x実測値の最小値と最大値の指定
  {
    if (getrange(xrange.first, xrange.second, opt[1]) == 0) // 範囲指定無
    { err = badopt(pname(), optstr, USAGE, -FDmax, FDmax); }
    else if (xrange.first >= xrange.second) // 範囲指定不適切
    { err = badopt(pname(), optstr, USAGE, -FDmax, FDmax); }
  }
  else if (opt[0] == "y")       // 目的変数の項目指定
  {
    char* ptr = nullptr;
    v_ys = (int)(strtol(opt[1].c_str(), &ptr, 10));
    if (v_ys < 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つ以上必要です！");
  int xe = v_xs + v_nt - 1, ye = v_ys + v_nt - 1;
  if (xe > (int)ni || ye > (int)ni) return disperr(pname(), "xとyの項目指定不適！");
  if (set_ngp() < 0) return -2; // -ngpオプションのセット
  if (set_res() < 0) return -3; // -resオプションのセット
  return 0;                     // エラーなし
}

int option::set_ngp()           // -ngpオプションのセット
{
  ngpmap["ce"] = ngpmap["data"] = ngpmap["ddata"] = ngpmap["line"] = ngpmap["lline"] = true;
  if (ngplst.size() > 0)
  {
    if (v_ngpname == "") return disperr(pname(), "-ngpnameオプションの指定が必要！");
    ngpmap["ce"] = ngpmap["data"] = ngpmap["ddata"] = ngpmap["line"] = ngpmap["lline"] = false;
    for (const string& ngpstr : ngplst)
    {
      mapsb_citr ngppos = ngpmap.find(ngpstr);
      if (ngppos == ngpmap.end())
      { return disperr(pname(), "-ngp=" + ngpstr + "オプション指定が不適切！"); }
      else
      { ngpmap[ngpstr] = true; }
    }
  }
  return 0;                     // エラーなし
}

int option::set_res()           // -resオプションのセット
{
  resmap["res"] = resmap["def"] = resmap["l"] = resmap["r"] = resmap["stat"] = true;
  if (reslst.size() > 0)
  {
    resmap["def"] = resmap["l"] = resmap["r"] = resmap["stat"] = false;
    for (const string& resstr : reslst)
    {
      mapsb_citr respos = resmap.find(resstr);
      if (resstr == "-")        // 全て出力せず
      { resmap["res"] = false; }
      else if (resstr[0] == '-' && resstr.size() > 1) // 特定の項目を出力せず
      {
        respos = resmap.find(resstr.substr(1));
        if (respos == resmap.end())
        { return disperr(pname(), "-res=" + resstr + "オプション指定が不適切！"); }
        else
        { resmap[resstr.substr(1)] = false; }
      }
      else if (resstr == "def" || resstr == "all")
      {
        resmap["def"] = resmap["l"] = resmap["r"] = resmap["stat"] = true;
      }
      else if (respos == resmap.end())
      { return disperr(pname(), "-res=" + resstr + "オプション指定が不適切！"); }
      else
      { resmap[resstr] = true; }
    }
  }
  return 0;                     // エラーなし
}

class linkcor                   // 連動係数用統計量クラス
{
private:
  int v_nd;                     // 差の例数
  double sigmadx;               // Σdx
  double sigmadx2;              // Σdx2
  double sigmady;               // Σdy
  double sigmady2;              // Σdy2
  double sigmadxy;              // Σdxdy
public:
  linkcor() : v_nd(0), sigmadx(0.0), sigmadx2(0.0), sigmady(0.0), sigmady2(0.0), sigmadxy(0.0)
  // デフォールトコンストラクタ
  {}
  linkcor(int nd0, double sigmadx0, double sigmadx20, double sigmady0, double sigmady20, double sigmadxy0) : v_nd(nd0), sigmadx(sigmadx0), sigmadx2(sigmadx20), sigmady(sigmady0), sigmady2(sigmady20), sigmadxy(sigmadxy0)
  // コンストラクタ
  {}
  linkcor(const linkcor& linkcor0) : v_nd(linkcor0.v_nd), sigmadx(linkcor0.sigmadx), sigmadx2(linkcor0.sigmadx2), sigmady(linkcor0.sigmady), sigmady2(linkcor0.sigmady2), sigmadxy(linkcor0.sigmadxy)
  // コピーコンストラクタ
  {}
  virtual ~linkcor()            // デストラクタ
  {}
  const linkcor& operator=(const linkcor& linkcor0); // コピー代入演算子
  const linkcor& operator+=(const linkcor& linkcor0); // 加算代入演算子
  int sigmad(double dx, double dy); // dxとdyの累積集計
  int nd() const                    // 差の例数を返す
  { return v_nd; }
  double mdx() const            // dxの平均を返す
  { return (nd() > 0) ? sigmadx / nd() : 0.0; }
  double mdy() const            // dyの平均を返す
  { return (nd() > 0) ? sigmady / nd() : 0.0; }
  double sdxx() const           // dxの単純平方和を返す
  { return sigmadx2; }
  double sdyy() const           // dyの単純平方和を返す
  { return sigmady2; }
  double sdxy() const           // dxとdyの単純積和を返す
  { return sigmadxy; }
  double sddx() const           // dxの標準変化を返す
  { return (nd() > 0) ? sqrt(sigmadx2 / nd()) : 0.0; }
  double sddy() const           // dyの標準変化を返す
  { return (nd() > 0) ? sqrt(sigmady2 / nd()) : 0.0; }
  double lamda() const;         // dxとdyの連動係数を返す
};

const linkcor& linkcor::operator=(const linkcor& linkcor0) // コピー代入演算子
{
  if (&linkcor0 == this) return *this; // 自分自身への代入を防ぐ
  v_nd = linkcor0.v_nd;  sigmadx = linkcor0.sigmadx;
  sigmadx2 = linkcor0.sigmadx2;  sigmady = linkcor0.sigmady;
  sigmady2 = linkcor0.sigmady2;  sigmadxy = linkcor0.sigmadxy;
  return *this;
}

const linkcor& linkcor::operator+=(const linkcor& linkcor0) // 加算代入演算子
{
  if (linkcor0.v_nd < 1) return *this;
  v_nd += linkcor0.v_nd;  sigmadx += linkcor0.sigmadx;
  sigmadx2 += linkcor0.sigmadx2;  sigmady += linkcor0.sigmady;
  sigmady2 += linkcor0.sigmady2;  sigmadxy += linkcor0.sigmadxy;
  return *this;
}

int linkcor::sigmad(double dx, double dy) // dxとdyの累積集計
{
  if (dx == BLANKf || dy == BLANKf) return 0; // 累積せず
  ++v_nd;  sigmadx += dx;  sigmadx2 += dx * dx;
  sigmady += dy;  sigmady2 += dy * dy;  sigmadxy += dx * dy;
  return 1;                     // 累積実行
}

double linkcor::lamda() const   // dxとdyの連動係数を返す
{
  double ssxy = sigmadx2 * sigmady2;
  return (ssxy > 0.0) ? sigmadxy / sqrt(ssxy) : 0.0;
}

double putlci(ostream& fout, double cc, int n, double l)
// 連動係数の信頼区間を出力する
{
  if (cc < 0.0 || 1.0 <= cc || n < 3 || fabs(l) == 1.0) return 0.0; // 計算できず
  double zl = log((1.0 + l) / (1.0 - l)) * 0.5;
  double ci = ptoz(1.0 - cc) / sqrt(double(n - 2));
  fout << "  連動係数の" << cc * 100.0 << "%信頼区間=tanh(" << zl << "±";
  fout << ci << ")(" << tanh(zl - ci) << "〜" << tanh(zl + ci) << ")\n";
  return ci;                    // 信頼幅を返す
}

class varclass                  // 変数クラス
{
private:
  int v_nc;                       // ケース数
  int v_nt;                       // 時期数
  pairdbl xrange;                 // xの最小値と最大値
  pairdbl dxrange;                // x変化量の最小値と最大値
  vector<statv<double>> xstatvec; // Xの全データ用基礎統計量ベクトル
  vector<statv<double>> ystatvec; // Yの全データ用基礎統計量ベクトル
  vector<statcor<double>> statcorvec; // 相関基礎統計量ベクトル
  vector<linkcor> linkcorvec;         // 連動係数用統計量ベクトル
public:
  varclass(int nt0) : v_nc(0), v_nt(nt0), xrange(FDmax, -FDmax), dxrange(FDmax, -FDmax), xstatvec(nt0 + 1), ystatvec(nt0 + 1), statcorvec(nt0 + 1), linkcorvec(nt0 + 1)
  // デフォールトコンストラクタ
  {}
  ~varclass()                   // デストラクタ
  {}
  int cal(option& opt, ldata& datatlst);
  // データを読み込み変数をセットする
  double xmin() const           // xの最小値を返す
  { return xrange.first; }
  double xmax() const           // xの最大値を返す
  { return xrange.second; }
  double dxmin() const          // x変化量の最小値を返す
  { return dxrange.first; }
  double dxmax() const          // x変化量の最大値を返す
  { return dxrange.second; }
  int nc() const                // ケース数を返す
  { return v_nc; }
  int nt() const                // 時期数を返す
  { return v_nt; }
  int xn(size_t idx) const      // xの例数を返す
  { return (idx < xstatvec.size()) ? xstatvec[idx].n() : 0; }
  int xnt() const               // xの全例数を返す
  { return xstatvec[nt()].n(); }
  double xmean(size_t idx) const // xの平均値を返す
  { return (idx < xstatvec.size()) ? xstatvec[idx].mean() : 0.0; }
  double xvr(size_t idx) const  // xの不偏分散を返す
  { return (idx < xstatvec.size()) ? xstatvec[idx].vr() : 0.0; }
  int yn(size_t idx) const      // yの例数を返す
  { return (idx < ystatvec.size()) ? ystatvec[idx].n() : 0; }
  int ynt() const               // yの全例数を返す
  { return ystatvec[nt()].n(); }
  double ymean(size_t idx) const // yの平均値を返す
  { return (idx < ystatvec.size()) ? ystatvec[idx].mean() : 0.0; }
  double yvr(size_t idx) const  // yの不偏分散を返す
  { return (idx < ystatvec.size()) ? ystatvec[idx].vr() : 0.0; }
  int cn(size_t idx) const      // xyペアの例数を返す
  { return (idx < statcorvec.size()) ? statcorvec[idx].n() : 0; }
  double cmx(size_t idx) const  // xyペアのxの平均値を返す
  { return (idx < statcorvec.size()) ? statcorvec[idx].mx() : 0.0; }
  double cmy(size_t idx) const  // xyペアのyの平均値を返す
  { return (idx < statcorvec.size()) ? statcorvec[idx].my() : 0.0; }
  double csxx(size_t idx) const // xyペアのxの平方和を返す
  { return (idx < statcorvec.size()) ? statcorvec[idx].sxx() : 0.0; }
  double csyy(size_t idx) const // xyペアのyの平方和を返す
  { return (idx < statcorvec.size()) ? statcorvec[idx].syy() : 0.0; }
  double csxy(size_t idx) const // xyペアの積和を返す
  { return (idx < statcorvec.size()) ? statcorvec[idx].sxy() : 0.0; }
  double cvxx(size_t idx) const // xyペアのxの不偏分散を返す
  { return (idx < statcorvec.size()) ? statcorvec[idx].vxx() : 0.0; }
  double cvyy(size_t idx) const // xyペアのyの不偏分散を返す
  { return (idx < statcorvec.size()) ? statcorvec[idx].vyy() : 0.0; }
  double cvxy(size_t idx) const // xyペアの共分散を返す
  { return (idx < statcorvec.size()) ? statcorvec[idx].vxy() : 0.0; }
  double crxy(size_t idx) const // xyペアの相関係数を返す
  { return (idx < statcorvec.size()) ? statcorvec[idx].rxy() : 0.0; }
  int lnd(size_t idx) const     // 変化量の例数を返す
  { return (idx < linkcorvec.size()) ? linkcorvec[idx].nd() : 0; }
  int lndt() const              // 変化量の全例数を返す
  { return linkcorvec[nt()].nd(); }
  double lmdx(size_t idx) const // x変化量の平均値を返す
  { return (idx < linkcorvec.size()) ? linkcorvec[idx].mdx() : 0.0; }
  double lmdy(size_t idx) const // y変化量の平均値を返す
  { return (idx < linkcorvec.size()) ? linkcorvec[idx].mdy() : 0.0; }
  double lsddx(size_t idx) const // x変化量の標準偏差を返す
  { return (idx < linkcorvec.size()) ? linkcorvec[idx].sddx() : 0.0; }
  double lsddy(size_t idx) const // y変化量の標準偏差を返す
  { return (idx < linkcorvec.size()) ? linkcorvec[idx].sddy() : 0.0; }
  double lsdxx(size_t idx) const // x変化量の平方和を返す
  { return (idx < linkcorvec.size()) ? linkcorvec[idx].sdxx() : 0.0; }
  double lsdyy(size_t idx) const // y変化量の平方和を返す
  { return (idx < linkcorvec.size()) ? linkcorvec[idx].sdyy() : 0.0; }
  double lsdxy(size_t idx) const // xy変化量の積和を返す
  { return (idx < linkcorvec.size()) ? linkcorvec[idx].sdxy() : 0.0; }
  double llamda(size_t idx) const // xy変化量の連動係数を返す
  { return (idx < linkcorvec.size()) ? linkcorvec[idx].lamda() : 0.0; }
};

int varclass::cal(option& opt, ldata& datalst)
// データを読み込み変数をセットする
{
  v_nc = 0;
  for (auto& cdatai : datalst)
  {
    double xb = cdatai.fdata(opt.xs() - 1); // Xの前値
    double yb = cdatai.fdata(opt.ys() - 1); // Yの前値
    bool bok = (xb != BLANKf && yb != BLANKf) ? true : false; // 前値の有無
    bool ncok = false;
    for (int i = 0; i < nt(); ++i)
    {
      double x = cdatai.fdata(opt.xs() + i - 1);
      double y = cdatai.fdata(opt.ys() + i - 1);
      if (x != BLANKf && y != BLANKf)
      {
        ncok = true;
        xstatvec[i].sigma(x);  ystatvec[i].sigma(y);
        statcorvec[i].sigma(x, y);
        if (x < xrange.first) xrange.first = x;
        if (x > xrange.second) xrange.second = x;
        if (i > 0 && bok)
        {
          double dx = x - xb, dy = y - yb;
          linkcorvec[i].sigmad(dx, dy);
          if (dx < dxrange.first) dxrange.first = dx;
          if (dx > dxrange.second) dxrange.second = dx;
        }
      }
    }
    if (ncok) ++v_nc;
  }
  if (v_nc < 1) return v_nc;    // 有効なケースなし
  for (int i = 0; i < nt(); ++i)
  {
    if (xstatvec[i].n() > 0)
    {
      xstatvec[nt()] += xstatvec[i];
      ystatvec[nt()] += ystatvec[i];
      statcorvec[nt()] += statcorvec[i];
      linkcorvec[nt()] += linkcorvec[i];
    }
  }
  return v_nc;                  // 有効なケース数を返す
}

class tcoregclass               // tcoregメインクラス
{
private:
  string fname;                 // ファイル名
  idd id;                       // データ識別情報
  ldata datalst;                // データリスト
  void put_resid(ostream& fout); // 分析結果のタイトルとデータ識別情報を出力する
  void put_resttid(ostream& fout, option& opt);
  // 時期合計分析結果のタイトルとデータ識別情報を出力する
  void put_coreg(ostream& fout, option& opt, int idx, const varclass& var);
  // 相関係数と回帰直線を出力する
  void put_link(ostream& fout, option& opt, int idx, const varclass& var);
  // 連動係数と連動回帰直線を出力する
  void plot_ngraph(option& opt, const varclass& var); // Ngraph用データ出力
  void plot_ngpid(ofstream& fout);
  // Ngraph用データのタイトルとデータ識別情報を出力する
  void plot_ngpttid(ofstream& fout, option& opt);
  // 時期合計Ngraph用データのタイトルとデータ識別情報を出力する
  void plot_ngpdata(ofstream& fout, option& opt, int idx);
  // Ngraph用に実測値のデータを出力する
  void plot_ngpline(ofstream& fout, option& opt, const varclass& var, int idx) const;
  // Ngraph用に回帰直線を出力する
  double set_xrange(double& xmin, double& xmax, option& opt, const varclass& var) const;
  // x実測値の最小値と最大値を設定して範囲を返す
  void plot_ngpce(ofstream& fout, option& opt, const varclass& var, int idx) const;
  // Ngraph用に信頼楕円を出力する
  void plot_ngpellipse(ofstream& fout, double mx, double my, double a, double b, double h) const;
  // 楕円のデータを出力する
  void plot_ngpddata(ofstream& fout, option& opt, int idx);
  // Ngraph用に変化量のデータを出力する
  void plot_ngplline(ofstream& fout, option& opt, const varclass& var, int idx) const;
  // Ngraph用に連動回帰直線を出力する
  double set_dxrange(double& dxmin, double& dxmax, option& opt, const varclass& var) const;
  // x変化量の最小値と最大値を設定して範囲を返す
public:
  tcoregclass() : fname(), id(), datalst() // コンストラクタ
  {}
  ~tcoregclass()                // デストラクタ
  {}
  int get(const char* c_fname, option& opt); // データを読み込み集計する
  void put(option& opt);                     // 基礎統計量を出力する
};

void tcoregclass::put_resid(ostream& fout)
// 分析結果のタイトルとデータ識別情報を出力する
{
  puttitle(fout, Title);
  fout << "\n";
  id.put_title(fout);
  id.put_stra(fout);
}

void tcoregclass::put_resttid(ostream& fout, option& opt)
// 時期合計分析結果のタイトルとデータ識別情報を出力する
{
  int inos = opt.xs(), inoe = opt.xs() + opt.nt() - 1;
  for (int i = inos; i <= inoe; ++i)
  {
    fout << "説明変数x" << setw(3) << i - opt.xs() + 1;
    fout << ":" << id.item(i - 1) << "\n";
  }
  inos = opt.ys();  inoe = opt.ys() + opt.nt() - 1;
  for (int i = inos; i <= inoe; ++i)
  {
    fout << "目的変数y" << setw(3) << i - opt.ys() + 1;
    fout << ":" << id.item(i - 1) << "\n";
  }
  fout << "\n";
}

void tcoregclass::put_coreg(ostream& fout, option& opt, int idx, const varclass& var)
// 相関係数と回帰直線を出力する
{
  if (opt.res("stat"))
  {
    fout << "x:";
    putmsd(fout, var.xn(idx), var.xmean(idx), var.xvr(idx));
    fout << "y:";
    putmsd(fout, var.yn(idx), var.ymean(idx), var.yvr(idx));
    putline(fout, LENline);
  }
  int n = var.cn(idx);
  int df = n - 2;
  double sxx = var.csxx(idx), syy = var.csyy(idx);
  double cxy = var.csxy(idx);
  double r = var.crxy(idx), b = (sxx > 0.0) ? cxy / sxx : 0.0;
  double a = var.cmy(idx) - var.cmx(idx) * b;
  double ci = 0.0;
  fout << "  相関係数:r=" << r << "  寄与率:r2乗=" << r * r << "\n";
  if (df > 0)
  {
    double zr = 0.0;
    ci = calrci(zr, opt.cc(), n, r);
    fout << "  相関係数の" << opt.cc() * 100.0 << "%信頼区間=tanh(" << zr << "±" << ci << ")(" << tanh(zr - ci) << "-" << tanh(zr + ci) << ")\n";
  }
  double vb = (df > 0 && sxx > 0.0) ? (syy - cxy * b) / (double)df / sxx : 0.0;
  fout << "  回帰直線:y=" << a;
  fout.setf(ios::showpos);
  fout << b << "x    回帰係数の標準誤差=";
  fout.unsetf(ios::showpos);
  fout << sqrt(vb) << "\n";
  if (df < 1) return;
  ci = calmcci(mcfis, opt.cc(), 1, 1, df, 1, vb);
  fout << "  回帰係数の" << opt.cc() * 100.0 << "%信頼区間=";
  fout << b << "±" << ci << "(" << b - ci << "-" << b + ci << ")\n";
  double vr = (df > 0 && fabs(r) < 1.0) ? (1.0 - r * r) / (double)df : 0.0;
  double t = (vr > 0.0) ? r / sqrt(vr) : 0.0;
  ldbl p = ttop(t, df);
  fout << "  相関係数・回帰係数の検定:t=" << t << "  自由度=" << df << "  有意確率:p=" << p << pind(p) << "\n";
}

void tcoregclass::put_link(ostream& fout, option& opt, int idx, const varclass& var)
// 連動係数と連動回帰直線を出力する
{
  if (opt.res("stat"))
  {
    fout << "xの変化量Δx:例数=" << left << setw(6) << var.lnd(idx);
    fout << "平均値=" << setw(12) << var.lmdx(idx);
    fout << "標準変化=" << var.lsddx(idx) << "\n";
    fout << "yの変化量Δy:例数=" << setw(6) << var.lnd(idx);
    fout << "平均値=" << setw(12) << var.lmdy(idx) << right;
    fout << "標準変化=" << var.lsddy(idx) << "\n";
    putline(fout, LENline);
  }
  double sx = var.lsdxx(idx), sy = var.lsdyy(idx), cxy = var.lsdxy(idx);
  double b = (sx > 0.0) ? cxy / sx : 0.0, l = var.llamda(idx);
  int df = var.lnd(idx) - 1;
  fout << "  連動係数:λ=" << l << "  寄与率:λ2乗=" << l * l << "\n";
  putlci(fout, opt.cc(), df + 1, l);
  double vb = (df > 0 && sx > 0.0) ? (sy - cxy * b) / (double)df / sx : 0.0;
  fout << "  連動回帰直線:Δy=" << b << "Δx    回帰係数の標準誤差=" << sqrt(vb) << "\n";
  if (df < 1) return;
  double ci = calmcci(mcfis, opt.cc(), 1, 1, df, 1, vb);
  fout << "  回帰係数の" << opt.cc() * 100.0 << "%信頼区間=";
  fout << b << "±" << ci << "(" << b - ci << "-" << b + ci << ")\n";
  double vl = (df > 0 && fabs(l) < 1.0) ? (1.0 - l * l) / double(df) : 0.0;
  double t = (vl > 0.0) ? l / sqrt(vl) : 0.0;
  ldbl p = ttop(t, df);
  fout << "  連動係数・回帰係数の検定:t=" << t << "  自由度=" << df << "  有意確率:p=" << p << pind(p) << "\n";
}

void tcoregclass::plot_ngraph(option& opt, const varclass& var)
// Ngraph用データ出力
{
  ofstream fout;
  if (ofopen(fout, (opt.ngpname()).c_str()) == false) return; // ファイルオープン失敗
  plot_ngpid(fout);
  int ntt = (opt.nt() > 1) ? opt.nt() : 0;
  bool delml = false;
  for (int i = 0; i < ntt; ++i)
  {
    if (var.xn(i) > 0)
    {
      if (opt.ti() || opt.ngp("data") || opt.ngp("ddata"))
      {
        if (delml) fout << NGPdelm << "\n";
        delml = false;
        fout << "%説明変数x:" << id.item(opt.xs() + i - 1) << "\n";
        fout << "%目的変数y:" << id.item(opt.ys() + i - 1) << "\n";
        if (opt.ngp("data"))    // 時期合計だけの時も出力する
        {
          plot_ngpdata(fout, opt, i);
          delml = true;
        }
        if (opt.ngp("line"))
        {
          if (delml) fout << NGPdelm << "\n";
          else delml = true;
          plot_ngpline(fout, opt, var, i);
        }
        if (opt.ngp("ce"))
        {
          if (delml) fout << NGPdelm << "\n";
          else delml = true;
          plot_ngpce(fout, opt, var, i);
        }
        if (i > 0)
        {
          if (opt.ngp("ddata")) // 時期合計だけの時も出力する
          {
            if (delml) fout << NGPdelm << "\n";
            else delml = true;
            plot_ngpddata(fout, opt, i);
          }
          if (opt.ngp("lline"))
          {
            if (delml) fout << NGPdelm << "\n";
            else delml = true;
            plot_ngplline(fout, opt, var, i);
          }
        }
      }
    }
  }
  if (var.xnt() < 1 || opt.tt() == false) return;
  if (opt.ngp("line") == false && opt.ngp("ce") == false && opt.ngp("lline") == false) return;
  if (delml) fout << NGPdelm << "\n";
  plot_ngpttid(fout, opt);
  delml = false;
  if (opt.ngp("line"))
  {
    plot_ngpline(fout, opt, var, ntt);
    delml = true;
  }
  if (opt.ngp("ce"))
  {
    if (delml) fout << NGPdelm << "\n";
    else delml = true;
    plot_ngpce(fout, opt, var, ntt);
  }
  if (opt.ngp("lline"))
  {
    if (delml) fout << NGPdelm << "\n";
    plot_ngplline(fout, opt, var, ntt);
  }
}

void tcoregclass::plot_ngpid(ofstream& fout)
// Ngraph用データのタイトルとデータ識別情報を出力する
{
  fout << "%";  puttitle(fout, Title);
  fout << "%\n";
  id.put_title(fout, formngp);
  id.put_stra(fout, formngp);
}

void tcoregclass::plot_ngpttid(ofstream& fout, option& opt)
// 時期合計Ngraph用データのタイトルとデータ識別情報を出力する
{
  int inos = opt.xs(), inoe = opt.xs() + opt.nt() - 1;
  for (int i = inos; i <= inoe; ++i)
  {
    fout << "%説明変数x" << setw(3) << i - opt.xs() + 1;
    fout << ":" << id.item(i - 1) << "\n";
  }
  inos = opt.ys();  inoe = opt.ys() + opt.nt() - 1;
  for (int i = inos; i <= inoe; ++i)
  {
    fout << "%目的変数y" << setw(3) << i - opt.ys() + 1;
    fout << ":" << id.item(i - 1) << "\n";
  }
  fout << "%時期合計\n";
}

void tcoregclass::plot_ngpdata(ofstream& fout, option& opt, int idx)
// Ngraph用に実測値のデータを出力する
{
  fout << "%実測値データ\n";
  fout << "%-s1 -dmark\n";
  for (auto& cdatai : datalst)
  {
    double x = cdatai.fdata(opt.xs() + idx - 1);
    double y = cdatai.fdata(opt.ys() + idx - 1);
    if (x != BLANKf && y != BLANKf) fout << x << " " << y << "\n";
  }
}

void tcoregclass::plot_ngpline(ofstream& fout, option& opt, const varclass& var, int idx) const
// Ngraph用に回帰直線を出力する
{
  fout << "%回帰直線\n";
  fout << "%-s1 -dline\n";
  double mx = var.cmx(idx), my = var.cmy(idx);
  double sxx = var.csxx(idx), syy = var.csyy(idx);
  double sxy = var.csxy(idx);
  double b = (sxx > 0.0) ? sxy / sxx : 0.0;
  double a = my - b * mx;
  double xmin = 0.0, xmax = 0.0;
  double xrange = set_xrange(xmin, xmax, opt, var);
  fout << xmin << " " << a + b * xmin << "\n";
  fout << xmax << " " << a + b * xmax << endl;
}

double tcoregclass::set_xrange(double& xmin, double& xmax, option& opt, const varclass& var) const
// x実測値の最小値と最大値を設定して範囲を返す
{
  if (opt.xmin() != BLANKf && opt.xmax() != BLANKf) // X軸の範囲指定あり
  { xmin = opt.xmin();  xmax = opt.xmax(); }
  else
  {
    xmin = var.xmin();  xmax = var.xmax();
    double xrange = xmax - xmin, xmg = 0.5;
    xmin -= xrange * xmg;  xmax += xrange * xmg;
  }
  return xmax - xmin;           // 範囲を返す
}

void tcoregclass::plot_ngpce(ofstream& fout, option& opt, const varclass& var, int idx) const
// Ngraph用に信頼楕円を出力する
{
  fout << "%-s1 -dpolygon\n";
  int n = var.cn(idx);
  double mx = var.cmx(idx), my = var.cmy(idx);
  double vxx = var.cvxx(idx), vyy = var.cvyy(idx), vxy = var.cvxy(idx);
  int dfr = n - 2;
  double cf = 1.0 / (double)n, rf = 1.0 + 1.0 / (double)n;
  double ccpr = opt.cc() * 100.0;
  fout << "%" << ccpr << "%信頼楕円\n";
  double a = 0.0, b = 0.0, h = 0.0;
  calelp(a, b, h, opt.cc(), 2, dfr, vxx * cf, vyy * cf, vxy * cf);
  plot_ngpellipse(fout, mx, my, a, b, h);
  fout << NGPdelm2 << "\n";
  fout << "%" << ccpr << "%許容楕円\n";
  calelp(a, b, h, opt.cc(), 2, dfr, vxx, vyy, vxy);
  plot_ngpellipse(fout, mx, my, a, b, h);
  fout << NGPdelm2 << "\n";
  fout << "%" << ccpr << "%棄却楕円\n";
  calelp(a, b, h, opt.cc(), 2, dfr, vxx * rf, vyy * rf, vxy * rf);
  plot_ngpellipse(fout, mx, my, a, b, h);
}

void tcoregclass::plot_ngpellipse(ofstream& fout, double mx, double my, double a, double b, double h) const
// 楕円のデータを出力する
{
  constexpr int RDdef = 60;     // 楕円の分割数
  constexpr double Dt = 2.0 * M_PI / (double)RDdef;
  double x0 = 0.0, y0 = 0.0, x = 0.0, y = 0.0, t = 0.0;
  for (int i = 0; i <= RDdef; ++i, t += Dt)
  {
    x0 = a * cos(t);  y0 = b * sin(t);
    x = mx + x0 * cos(h) - y0 * sin(h);  y = my + x0 * sin(h) + y0 * cos(h);
    fout << x << " " << y << "\n";
  }
}

void tcoregclass::plot_ngpddata(ofstream& fout, option& opt, int idx)
// Ngraph用に変化量のデータを出力する
{
  fout << "%変化量データ\n";
  fout << "%-s1 -dmark\n";
  for (auto& cdatai : datalst)
  {
    double xb = cdatai.fdata(opt.xs() - 1); // Xの前値
    double yb = cdatai.fdata(opt.ys() - 1); // Yの前値
    if (xb != BLANKf && yb != BLANKf)
    {
      double x = cdatai.fdata(opt.xs() + idx - 1);
      double y = cdatai.fdata(opt.ys() + idx - 1);
      if (x != BLANKf && y != BLANKf) fout << x - xb << " " << y - yb << "\n";
    }
  }
}

void tcoregclass::plot_ngplline(ofstream& fout, option& opt, const varclass& var, int idx) const
// Ngraph用に連動回帰直線を出力する
{
  fout << "%連動回帰直線\n";
  double sx = var.lsdxx(idx), sy = var.lsdyy(idx), cxy = var.lsdxy(idx);
  double b = (sx > 0.0) ? cxy / sx : 0.0;
  double dxmin = 0.0, dxmax = 0.0;
  double dxrange = set_dxrange(dxmin, dxmax, opt, var);
  fout << dxmin << " " << b * dxmin << "\n";
  fout << dxmax << " " << b * dxmax << "\n";
}

double tcoregclass::set_dxrange(double& dxmin, double& dxmax, option& opt, const varclass& var) const
// x変化量の最小値と最大値を設定して範囲を返す
{
  if (opt.dxmin() != BLANKf && opt.dxmax() != BLANKf) // X軸の範囲指定あり
  { dxmin = opt.dxmin();  dxmax = opt.dxmax(); }
  else
  {
    dxmin = var.dxmin();  dxmax = var.dxmax();
    double dxrange = dxmax - dxmin, dxmg = 0.5;
    dxmin -= dxrange * dxmg;  dxmax += dxrange * dxmg;
  }
  return dxmax - dxmin;         // 範囲を返す
}

int tcoregclass::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 tcoregclass::put(option& opt) // 集計結果を出力する
{
  varclass var(opt.nt());           // 計算用変数オブジェクト
  var.cal(opt, datalst);
  if (opt.res("res"))           // 計算結果出力
  {
    put_resid(cout);
    if (var.nc() < 1)
    {
      cout << "説明変数x:" << id.item(opt.xs() - 1) << "\n";
      cout << "目的変数y:" << id.item(opt.ys() - 1) << "\n";
      cout << "※有効なデータなし！\n\n" << endl;
      return;
    }
    int ntt = (opt.nt() > 1) ? opt.nt() : 0;
    bool putresi = true, putrest = false;
    for (int i = 0; i < ntt; ++i)
    {
      if (var.xn(i) > 0)
      {
        if (opt.ti())
        {
          cout << "説明変数x:" << id.item(opt.xs() + i - 1) << "\n";
          cout << "目的変数y:" << id.item(opt.ys() + i - 1) << "\n";
          putresi = false;
          if (opt.res("r"))
          {
            putline(cout, LENline);
            put_coreg(cout, opt, i, var);
            putresi = true;
          }
          if (opt.res("l") && var.lnd(i) > 0)
          {
            putline(cout, LENline);
            put_link(cout, opt, i, var);
            putresi = true;
          }
          if (putresi) putline(cout, LENline);
          cout << "\n";
        }
      }
    }
    if (var.xnt() > 0 && opt.tt())
    {
      if (opt.ti() == false) put_resttid(cout, opt);
      cout << "・時期合計\n";
      putrest = false;
      if (opt.res("r"))
      {
        putline(cout, LENline);
        put_coreg(cout, opt, ntt, var);
        putrest = true;
      }
      if (opt.res("l") && var.lndt() > 0)
      {
        putline(cout, LENline);
        put_link(cout, opt, ntt, var);
        putrest = true;
      }
    }
    if (putrest)
    {
      putline(cout, LENline);
      cout << "\n";
    }
    else if (putresi == false)
    { cout << "\n"; }
    cout << endl;
  }
  if (var.nc() > 0 && opt.ngpname() != "") plot_ngraph(opt, var); // Ngraph用データ出力
}

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