// clust.C Version 7.45  Copyleft (c) 杉本典夫
// Time-stamp: <2026-04-12 20:59:43 tac>
//
// クラスター分析を行う。
//
// 使用ライブラリ:libdans.so, libfun.so
// コンパイル法 :g++ -O -s -W clust clust.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 =            // プログラム使用法表示用
{
  "clust  Version 7.45          Copyleft(Vo\\oV) 2026 杉本典夫",
  "書式:clust  [<オプション> …]  [<データファイル> …]",
  "機能:最初の項目をケースID,以後の項目をデータとしたクラスター分析を行うのだ。",
  "<オプション>",
  "  -den[=<樹形図の深さ>]:樹形図の深さを指定するよ。(-den:-den=2)",
  "  -dis={euc|euc2|md|md2|seuc|seuc2}:距離の種類を指定するよ。",
  "    euc:エウクレイデス距離  euc2:エウクレイデス平方距離(オプション省略時)",
  "    md:マハラノビス汎距離  md2:マハラノビス平方距離",
  "    seuc:標準化エウクレイデス距離  seuc2:標準化エウクレイデス平方距離",
  "  -dnsname[=<出力ファイル名>]:DANSデータを出力するファイル名を指定するよ。",
  "    オプション省略時:DANSデータ出力せず  -dnsname:-dnsname=clust_dns.d",
  "    出力項目:ケースごとのIDと順位とクラスター番号",
  "  -f=<ファイル>:指定したファイルからオプションを読み込むよ。",
  "  -h[=<コマンド名>]:オンラインマニュアルを表示するよ。(-h:-h=man)",
  "  -meth={cent|fur|group|med|near|ward}:融合法を指定するよ。",
  "    cent:重心法  fur:最長距離法  group:群平均法  med:メジアン法",
  "    near:最短距離法  ward:ウォード法(オプション省略時)",
  "  -nc=[<最大クラスター数>][-[<最小クラスター数>]]:出力するクラスター数の指定。",
  "  -ngpname[=<出力ファイル名>]:Ngraph用樹形図の出力ファイル名を指定するよ。",
  "    オプション省略時:出力せず  -ngpname:-ngpname=clust_ngp.dat",
  "  -res=[-][{def|all|den|sil|stat|v}][,…]:結果の出力項目を指定するよ。",
  "    -:-を付けた項目を出力しない(-res=-:全ての項目を出力しない)",
  "    def:オプション省略時と同じ項目を出力する  all:全ての項目を出力する",
  "    den:樹形図  sil:シルエット分析  stat:基礎統計量  v:計算過程"
};

const string Title = "=== クラスター分析(cluster analysis) ==="; // 表題
const string DNSNAMEdef = "clust_dns.d"; // 省略時DANSファイル名
const string NGPNAMEdef = "clust_ngp.dat"; // 省略時Ngraphファイル名
constexpr int DENdef = 2;                  // 省略時樹形図の深さ

enum distype                    // 距離の種類
{
  disnon,                       // 指定なし
  diseuc,                       // エウクレイデス距離
  diseuc2,                      // エウクレイデス平方距離
  dismd,                        // マハラノビス汎距離
  dismd2,                       // マハラノビス平方距離
  disseuc,                      // 標準化エウクレイデス距離
  disseuc2,                     // 標準化エウクレイデス平方距離
  disetc                        // その他
};

enum methtype                   // 融合法の種類
{
  methnon,                      // 指定なし
  methcent,                     // 重心法
  methfur,                      // 最長距離法
  methgroup,                    // 群平均法
  methmed,                      // メジアン法
  methnear,                     // 最短距離法
  methward,                     // ウォード法
  methetc                       // その他
};

class option                    // オプションクラス
{
private:
  string v_path;                // プログラムのパス
  string v_pname;               // プログラム名
  distype v_dis;                // 距離の種類
  string v_disstr;              // 距離の種類(文字列型)
  methtype v_meth;              // 融合法の種類
  string v_methstr;             // 融合法の種類(文字列型)
  int v_den;                    // 樹形図の深さ
  string v_dnsname;             // DANSファイル名
  string v_ngpname;             // Ngraphファイル名
  pairint ncrange;              // 出力クラスター数
  lstr reslst;                  // 分析結果出力リスト
  mapsb resmap;                 // 分析結果出力の有無マップ
public:
  option(const char* c_pname = nullptr); // コンストラクタ
  ~option()                              // デストラクタ
  {}
  const string& path() const    // プログラムのパスを返す
  { return v_path; }
  const string& pname() const   // プログラム名を返す
  { return v_pname; }
  distype dis() const           // 距離の種類を返す
  { return v_dis; }
  const string& disstr();       // 距離の種類を文字列で返す
  int den() const               // 樹形図の深さを返す
  { return v_den; }
  const string& dnsname() const // DANSデータファイル名を返す
  { return v_dnsname; }
  methtype meth() const         // 融合法の種類を返す
  { return v_meth; }
  const string& methstr();      // 融合法の種類を文字列で返す
  const string& ngpname() const // Ngraphデータファイル名を返す
  { return v_ngpname; }
  int ncmax() const             // 出力する最大クラスター数を返す
  { return ncrange.first; }
  int ncmin() const             // 出力する最小クラスター数を返す
  { return ncrange.second; }
  bool res(const string& code); // 分析結果出力フラッグを返す
  void setdef()                 // 省略時オプション指定を設定する
  { getoptdef(*this); }
  void set(const string& optstr, char devc = OPTdev);
  // オプション指定によってオプションを設定する
  int check(idd& id);           // オプションのチェックとセット
  int set_res();                // -resオプションのセット
};

option::option(const char* c_pname) : v_path(), v_pname(), v_dis(diseuc2), v_disstr(), v_den(DENdef), v_dnsname(), v_meth(methward), v_methstr(), v_ngpname(), ncrange(INT_MAX, 1), reslst(), resmap()
// コンストラクタ
{
  if (c_pname == nullptr) return; // 初期設定せず
  v_pname = getname(c_pname);     // パス名なしのプログラム名取得
  v_path = c_pname;
  strsize_t idx = v_path.find(v_pname);
  if (idx != string::npos) v_path.erase(idx); // プログラム名削除
}

bool option::res(const string& code) // 分析結果出力フラッグを返す
{
  mapsb_citr resitr = resmap.find(code);
  if (resitr == resmap.end()) return false; // code無し！
  return resmap[code];                     // true/falseを返す
}

void option::set(const string& optstr, char devc)
// オプション指定によってオプションを設定する
{
  string opt[2];                // オプションのopt[0]:項目,opt[1]:内容
  if (getopt(opt, optstr, devc) < 1) return; // オプション指定なし
  int err = 0;
  if (opt[0] == "dis")          // 距離の種類指定
  {
    if (opt[1].empty() || opt[1] == "euc2") v_dis = diseuc2;
    else if (opt[1] == "euc") v_dis = diseuc;
    else if (opt[1] == "md") v_dis = dismd;
    else if (opt[1] == "md2") v_dis = dismd2;
    else if (opt[1] == "seuc") v_dis = disseuc;
    else if (opt[1] == "seuc2") v_dis = disseuc2;
    else err = badopt(pname(), optstr, USAGE);
  }
  else if (opt[0] == "den")     // 樹形図の深さ指定
  {
    if (opt[1].empty())
    { v_den = DENdef; }
    else
    {
      char* ptr = nullptr;
      v_den = strtol(opt[1].c_str(), &ptr, 0);
      if (v_den < 1) err = badopt(pname(), optstr, USAGE, 1.0, (double)INT_MAX);
    }
  }
  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] == "f")       // オプション指定ファイル
  { err = fgetopt(*this, opt[1], true); }
  else if (opt[0] == "h")       // オンラインヘルプ表示
  { err = help(pname(), USAGE, opt[1]); }
  else if (opt[0] == "meth")    // 融合法の種類指定
  {
    if (opt[1].empty() || opt[1] == "ward") v_meth = methward;
    else if (opt[1] == "cent") v_meth = methcent;
    else if (opt[1] == "fur") v_meth = methfur;
    else if (opt[1] == "group") v_meth = methgroup;
    else if (opt[1] == "med") v_meth = methmed;
    else if (opt[1] == "near") v_meth = methnear;
    else err = badopt(pname(), optstr, USAGE);
  }
  else if (opt[0] == "nc")      // 出力するクラスター数指定
  {
    if (getrange(ncrange.first, ncrange.second, opt[1]) < 1) err = badopt(pname(), optstr, USAGE, INT_MAX, 1);
    if (ncrange.first < 1 || ncrange.second < 1 || ncrange.first < ncrange.second) err = badopt(pname(), optstr, USAGE, INT_MAX, 1);
  }
  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                          // 間違ったオプション指定
  { err = nonopt(pname(), optstr, USAGE); }
  if (err > 0) throw help_except("ヘルプ表示");
  if (err < 0) throw option_error(string("オプション指定違反→") + optstr);
}

int option::check(idd& id)      // オプションのチェックとセット
{
  if (id.ni() < 2) return disperr(pname(), "項目が2個以上必要だよ！");
  if (set_res() < 0) return -3; // -resオプションのエラー
  return 0;                     // エラーなし
}

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

const string& option::disstr()  // 距離の種類を文字列で返す
{
  if (v_disstr.empty())
  {
    switch (dis())
    {
    case diseuc : v_disstr = "エウクレイデス距離";
      break;
    case diseuc2 : v_disstr = "エウクレイデス平方距離";
      break;
    case dismd : v_disstr = "マハラノビス汎距離";
      break;
    case dismd2 : v_disstr = "マハラノビス平方距離";
      break;
    case disseuc : v_disstr = "標準化エウクレイデス距離";
      break;
    case disseuc2 : v_disstr = "標準化エウクレイデス平方距離";
      break;
    default : v_disstr = "その他";
    }
  }
  return v_disstr;
}

const string& option::methstr() // 融合法の種類を文字列で返す
{
  if (v_methstr.empty())
  {
    switch (meth())
    {
    case methcent : v_methstr = "重心法(centroid method)";
      break;
    case methfur : v_methstr = "最長距離法(furthest neighbour method)";
      break;
    case methgroup : v_methstr = "群平均法(group-average method)";
      break;
    case methmed : v_methstr = "メジアン法(median method)";
      break;
    case methnear : v_methstr = "最短距離法(nearest neighbour method)";
      break;
    case methward : v_methstr = "ウォード法(Ward's method)";
      break;
    default : v_methstr = "その他";
    }
  }
  return v_methstr;
}

class cnode                     // クラスター用ノードクラス
{
private:
  string v_cid;                 // ケースID
  int v_lno;                    // 最初に融合された左クラスター番号
  int v_rno;                    // 最初に融合された右クラスター番号
  int v_lidx;                   // 最初に融合された左クラスターのインデックス
  int v_ridx;                   // 最初に融合された右クラスターのインデックス
public:
  cnode() : v_cid(), v_lno(-1), v_rno(-1), v_lidx(-1), v_ridx(-1)
  // デフォールトコンストラクタ
  {}
  virtual ~cnode()              // デストラクタ
  {}
  const cnode& operator=(const cnode& cnode0); // コピー代入演算子
  const string& cid() const                    // ケースIDを返す
  { return v_cid; }
  const string& cid(const string& cid0) // ケースIDを設定する
  { return (v_cid = cid0); }
  int lno() const               // 左クラスターの番号を返す
  { return v_lno; }
  int lno(int lno0)             // 左クラスターの番号を設定する
  { return (v_lno = lno0); }
  int rno() const               // 右クラスターの番号を返す
  { return v_rno; }
  int rno(int rno0)             // 右クラスターの番号を設定する
  { return (v_rno = rno0); }
  int lidx() const              // 左クラスターのインデックスを返す
  { return v_lidx; }
  int lidx(int lidx0)          // 左クラスターのインデックスを設定する
  { return (v_lidx = lidx0); }
  int ridx() const              // 右クラスターのインデックスを返す
  { return v_ridx; }
  int ridx(int ridx0)          // 右クラスターのインデックスを設定する
  { return (v_ridx = ridx0); }
};

const cnode& cnode::operator=(const cnode& cnode0) // コピー代入演算子
{
  if (&cnode0 == this) return *this; // 自分自身への代入を防ぐ
  v_cid = cnode0.v_cid;  v_lno = cnode0.v_lno;  v_rno = cnode0.v_rno;
  v_lidx = cnode0.v_lidx;  v_ridx = cnode0.v_ridx;
  return *this;
}

class varclass                 // 変数クラス
{
private:
  vector<statv<double>> statvec; // 基礎統計量ベクトル
  mat<double> v_smat;            // 積和行列
  mat<double> gmat;              // 計量行列
  vector<cnode> ccvec;           // クラスター用ノードベクトル
  mat<double> dcmat0;            // クラスター間の初期距離行列
  mat<double> dcmat;             // クラスター間の距離行列
  vecdbl v_mevec;                // 平均値ベクトル
  vecdbl v_vrvec;                // 分散ベクトル
  int v_ni;                      // 項目数
  restype v_resvar;              // 処理結果の種類
  void set_cnode(option& opt, ldata& datalst);
  // クラスター分析用ノードを設定する
  double cal_dis(option& opt, datac& cdatai, datac& cdataj) const;
  // 個体間距離を返す
  int set_ctree(vecint& cnovec, int rootidx, int cidx);
  // クラスターの個体番号を再帰的にセットして個体数を返す
public:
  varclass(size_t ni) : statvec(ni - 1), v_smat(ni, ni), gmat(ni - 1, ni - 1), ccvec(), dcmat0(), dcmat(), v_mevec(ni, 0.0), v_vrvec(ni, 0.0), v_ni(ni), v_resvar(resnon)
  // コンストラクタ
  {}
  ~varclass()                   // デストラクタ
  {}
  cnode& cc(int idx)            // クラスター用ノードを返す
  { return (0 <= idx && idx < nc() * 2 - 1) ? ccvec[idx] : ccvec[0]; }
  double dc0(int i, int j) const // クラスター間の初期距離を返す
  { return (0 <= i && i < nc() && 0 <= j && j < nc()) ? dcmat0(i, j) : 0.0; }
  double dc(int i, int j) const // クラスター間の距離を返す
  { return (0 <= i && i < nc() && 0 <= j && j < nc()) ? dcmat(i, j) : 0.0; }
  const mat<double>& smat() const // 積和行列を返す
  { return v_smat; }
  const vecdbl& mevec() const   // 平均値ベクトルを返す
  { return v_mevec; }
  const vecdbl& vrvec() const   // 分散ベクトルを返す
  { return v_vrvec; }
  int ni() const                // 項目数を返す
  { return v_ni; }
  int np() const                // パラメータ数を返す
  { return ni() - 1; }
  int nc() const                // ケース数を返す
  { return statvec[0].n(); }
  restype resvar() const        // 処理結果の種類を返す
  { return v_resvar; }
  void cal(option& opt, ldata& datalst); // データを集計して各種統計量を求める
  void fusion_clust(int& g, int& h, int& ng, int& nh, int rootidx);
  // クラスターを融合する
  void update_clust(option& opt, int g, int h, int ng, int nh);
  // クラスター間距離を更新する
  int mesil(ostream& fout, option& opt, double& mes);
  // シルエット係数平均値を計算する
  void put_dcmat(ostream& fout);
  // dcmatの内容を出力する
  void put_ccvec(ostream& fout) const;
  // ccvecの内容を出力する
};

void varclass::cal(option& opt, ldata& datalst)
// データを集計して各種統計量を求める
{
  int dim = np() + 1;
  vecdbl fdatavec(dim, 0.0);
  fdatavec[0] = 1.0;
  for (int j = 0; j < np(); ++j) statvec[j].clear();
  for (auto& cdatai : datalst)
  {
    bool blank = false;
    for (int j = 1; blank == false && j < ni(); ++j)
    {
      fdatavec[j] = cdatai.fdata(j);
      if (fdatavec[j] == BLANKf) blank = true;
    }
    if (cdatai.sdata(0) == BLANKs || blank)
    { cdatai.assign(0, BLANKs); } // 以後のチェックのためにIDをブランクにする
    else
    {
      for (int i = 0; i < dim; ++i)
      {
        if (i < np()) statvec[i].sigma(fdatavec[i + 1]);
        double xi = fdatavec[i];
        for (int j = i; j < dim; ++j) v_smat.add(i, j, xi * fdatavec[j]);
      }
    }
  }
  int nc = statvec[0].n();
  if (nc < 2)
  {
    v_resvar = 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);              // 積和行列にする
  switch (opt.dis())
  {
  case dismd :
  case dismd2 : gmat = v_smat.sub(1, 1, np(), np());
    gmat *= 1.0 / (double)(nc - 1); // 共分散行列にする
    gmat.inv();                     // 逆行列にする
    break;
  case disseuc :
  case disseuc2 :
    for (int i = 0; i < np(); ++i) gmat.assign(i, i, 1.0 / statvec[i].vr());
    break;
  default : for (int i = 0; i < np(); ++i) gmat.assign(i, i, 1.0);
  }
  dcmat0.resize(nc, nc);  dcmat.resize(nc, nc);  ccvec.resize(nc * 2 - 1);
  v_mevec[0] = (double)nc;  v_vrvec[0] = (double)(nc - 1);
  for (int i = 1; i < ni(); ++i)
  { v_mevec[i] = statvec[i - 1].mean();  v_vrvec[i] = statvec[i - 1].vr(); }
  set_cnode(opt, datalst);
  v_resvar = resok;
}

void varclass::set_cnode(option& opt, ldata& datalst)
// クラスター分析用ノードを設定する
{
  size_t i = 0;
  ldata_itr ditri = datalst.begin();
  while (i < dcmat.r())
  {
    if (ditri->sdata(0) != BLANKs)
    {
      ccvec[i].cid(ditri->sdata(0));
      ccvec[i].ridx(i);
      dcmat0.assign(i, i, 0.0);
      dcmat.assign(i, i, 1.0);
      ldata_itr ditrj = ditri;
      ++ditrj;
      size_t j = i + 1;
      while (j < dcmat.c())
      {
        if (ditrj->sdata(0) != BLANKs)
        {
          dcmat.assign(i, j, cal_dis(opt, *ditri, *ditrj));
          dcmat.assign(j, i, dcmat(i, j));
          dcmat0.assign(i, j, dcmat(i, j));
          dcmat0.assign(j, i, dcmat(i, j));
          ++j;
        }
        ++ditrj;
      }
      ++i;
    }
    ++ditri;
  }
}

double varclass::cal_dis(option& opt, datac& cdatai, datac& cdataj) const
// 個体間距離を返す
{
  vec<double> xv(np());
  for (int i = 0; i < np(); ++i) xv.assign(i, cdatai.fdata(i + 1) - cdataj.fdata(i + 1));
  vec<double> tv(xv.trn() * gmat);
  double dis = inp(tv, xv);
  if (opt.dis() == diseuc2 || opt.dis() == dismd2 || opt.dis() == disseuc2) return dis;
  return sqrt(dis);
}

void varclass::fusion_clust(int& g, int& h, int& ng, int& nh, int rootidx)
// クラスターを融合する
{
  double dmin = FDmax;
  for (size_t k = 0; k < dcmat.r(); ++k)
  {
    if (dcmat(k, k) > 0.0)     // 融合されていないクラスター
    {
      int imin = k;
      double mink = FDmax;
      for (size_t j = k + 1; j < dcmat.c(); ++j)
      {
        if (dcmat(j, j) > 0.0)
        {
          if (dcmat(k, j) < mink) { mink = dcmat(k, j);  imin = j; }
        }
      }
      if (mink < dmin) { dmin = mink;  g = k;  h = imin; }
    }
  }
  ng = (int)(dcmat(g, g));  nh = (int)(dcmat(h, h));
  dcmat.assign(g, g, (double)(ng + nh));  dcmat.assign(h, h, 0.0);
  ccvec[rootidx].cid(ccvec[g].cid());
  ccvec[rootidx].lno(g);  ccvec[rootidx].rno(h);
  ccvec[rootidx].lidx(ccvec[g].ridx());
  ccvec[rootidx].ridx(ccvec[h].ridx());
  ccvec[g].ridx(rootidx);
}

void varclass::update_clust(option& opt, int g, int h, int ng, int nh)
// クラスター間距離を更新する
{
  for (size_t k = 0; k < dcmat.r(); ++k)
  {
    if (dcmat(k, k) > 0.0 && (int)k != g)
    {
      double dgk = dcmat(g, k), dhk = dcmat(h, k), dgh = dcmat(g, h);
      double dis = 0.0, ngh = (double)(ng + nh), nk = dcmat(k, k);
      switch (opt.meth())
      {
      case methnear : dis = (dgk < dhk ? dgk : dhk);
        break;
      case methfur : dis = (dgk > dhk ? dgk : dhk);
        break;
      case methmed : dis = (dgk + dhk) * 0.5 - dgh * 0.25;
        break;
      case methcent : dis = dgk * dgk * ng + dhk * dhk * nh - dgh * dgh * ng * nh / ngh;
        dis = sqrt(dis / ngh);
        break;
      case methgroup : dis = sqrt(dgk * dgk / ngh * ng + dhk * dhk / ngh * nh);
        break;
      case methward : dis = dgk * dgk * (ng + nk) + dhk * dhk * (nh + nk) - dgh * dgh * nk;
        dis = sqrt(dis / (ngh + nk));
        break;
      default : dis = 0.0;
      }
      dcmat.assign(k, g, dis);  dcmat.assign(g, k, dis);
    }
  }
}

int varclass::mesil(ostream& fout, option& opt, double& mes)
// シルエット係数平均値を計算する
{
  mes = 0.0;
  int nclust = 0, nc2 = 0;
  vecint civec(nc(), 0), cjvec(nc(), 0);
  for (int i = 0; i < nc(); ++i)
  {
    if (dc(i, i) > 0.0)         // 融合されていないクラスター
    {
      ++nclust;
      int ni = set_ctree(civec, cc(i).ridx(), 0);
      if (ni > 1)               // 個体が2つ以上のクラスター
      {
        for (int k = 0; k < ni; ++k)
        {
          int ck = civec[k];
          double ai = 0.0;
          for (int j = 0; j < ni; ++j)
          {
            if (j != k) ai += dc0(ck, civec[j]);
          }
          ai /= (double)(ni - 1); // クラスター内平均距離
          double dismin = FDmax;
          int idxmin = ck;
          for (int j = 0; j < nc(); ++j)
          {
            if (dc(j, j) > 0.0 && j != i && j != ck)
            {
              if (dc0(ck, j) < dismin) { dismin = dc0(ck, j);  idxmin = j; }
            }
          }
          if (idxmin != ck)
          {
            int nj = set_ctree(cjvec, cc(idxmin).ridx(), 0);
            double bi = 0.0;
            for (int j = 0; j < nj; ++j) bi += dc0(ck, cjvec[j]);
            bi /= (double)nj;   // 最近隣クラスター間平均距離
            double sil = (bi - ai) / max(ai, bi);
            mes += sil;  ++nc2;
            if (opt.res("v"))
            {
              fout << "  " << "クラスター" << cc(i).cid() << ":" << cc(ck).cid() << "-最近隣:" << cc(idxmin).cid();
              fout << " sil=" << setw(10) << sil << " mes=" << setw(10) << mes / (double)(nc()) << "\n";
            }
          }
        }
      }
    }
  }
  mes = mes / (double)(nc());
  return nclust;                // クラスター数を返す
}

int varclass::set_ctree(vecint& cnovec, int rootidx, int cidx)
// クラスターの個体番号を再帰的にセットして個体数を返す
{
  if (cc(rootidx).lidx() < 0)
  {
    cnovec[cidx] = rootidx;
    return cidx + 1;
  }
  cidx = set_ctree(cnovec, cc(rootidx).lidx(), cidx);
  cidx = set_ctree(cnovec, cc(rootidx).ridx(), cidx);
  return cidx;
}

void varclass::put_dcmat(ostream& fout)
// dcmatの内容を出力する
{
  fout << setw(6) << " ";
  for (size_t j = 0; j < dcmat.c(); ++j) fout << setw(9) << cc(j).cid();
  fout << "\n";
  for (size_t i = 0; i < dcmat.r(); ++i)
  {
    fout << setw(6) << cc(i).cid();
    for (size_t j = 0; j < dcmat.c(); ++j) fout << setw(9) << dcmat(i, j);
    fout << "\n";
  }
}

void varclass::put_ccvec(ostream& fout) const
// ccvecの内容を出力する
{
  for (size_t i = 0; i < ccvec.size(); ++i)
  {
    fout << "  ccvec[" << setw(2) << i << "]=" << setw(3) << ccvec[i].cid();
    fout << ":lno=" << setw(2) << ccvec[i].lno() << " rno=" << setw(2) << ccvec[i].rno();
    fout << " lidx=" << setw(2) << ccvec[i].lidx() << " ridx=" << setw(2) << ccvec[i].ridx() << "\n";
  }
}

class clustclass              // clustメインクラス
{
private:
  string fname;                 // ファイル名
  idd id;                       // データ識別情報
  ldata datalst;                // データリスト
  void put_resid(ostream& fout);
  // 解析結果のタイトルとデータ識別情報を出力する
  int put_clust(ostream& fout, option& opt, varclass& var);
  // クラスター分析の結果を出力する
  int put_ctree(ostream& fout, int rootidx, int jidx, varclass& var);
  // クラスターのケースIDを再帰的に出力する
  void put_cid(ostream& fout, option& opt, varclass& var);
  // クラスターのIDを出力する
  void put_dns(option& opt, int rootidx, varclass& var);
  // ケースIDとクラスター番号をデータファイル形式で出力する
  int put_dnsctree(ofstream& fout, datac& cdata, int& gno, int tlevel0, int tlevel, int rootidx, int jidx, varclass& var);
  // ケースIDとクラスター番号を再帰的に出力する
  void plot_dend(option& opt, int rootidx, varclass& var);
  // 樹形図用データを出力する
  void plot_dtree(ofstream& fout, double& xc, double& yc, int rootidx, varclass& var) const;
  // 樹形図の幹データを再帰的に出力する
public:
  clustclass() : fname(), id(), datalst() // コンストラクタ
  {}
  ~clustclass()                 // デストラクタ
  {}
  int get(const char* c_fname, option& opt); // データを読み込み累積する
  void put(option& opt);                     // 集計結果を出力する
};

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

int clustclass::put_clust(ostream& fout, option& opt, varclass& var)
// クラスター分析の結果を出力する
{
  if (opt.res("v"))
  {
    fout << "・ステップ0:ccvec初期値とdcmat初期値\n";
    var.put_ccvec(fout);
    var.put_dcmat(fout);
    fout << "  シルエット係数平均値=0(クラスター数=" << var.nc() << ")\n\n";
  }
  ldata_itr ditr = datalst.begin();
  int i = 0, rootidx = 0;
  while (i < var.nc() - 1)
  {
    if (ditr->sdata(0) != BLANKs)
    {
      int g = 0, h = 0, ng = 0, nh = 0;
      rootidx = var.nc() + i;
      var.fusion_clust(g, h, ng, nh, rootidx);
      int nc = var.nc() - i - 1;
      bool putd = false;
      bool putnc = (opt.ncmin() <= nc && nc <= opt.ncmax()) ? true : false;
      if (opt.res("den") && putnc)
      {
        fout << "・ステップ" << i + 1 << ":[" << (var.cc(g)).cid() << "]←[" << (var.cc(h)).cid() << "]融合  クラスター間の距離=" << var.dc(g, h) << "\n";
        put_cid(fout, opt, var);
        putd = true;
      }
      var.update_clust(opt, g, h, ng, nh);
      if (opt.res("v") && putnc)
      {
        var.put_ccvec(fout);
        var.put_dcmat(fout);
        putd = true;
      }
      if (opt.res("sil") && putnc)
      {
        double mes = 0.0;
        int nclust = var.mesil(fout, opt, mes);
        fout << "  シルエット係数平均値=" << mes << "(クラスター数=" << nclust << ")\n";
        putd = true;
      }
      if (putd) { fout << "\n"; }
      ++i;
    }
    ++ditr;
  }
  return rootidx;               // ルートインデックスを返す
}

int clustclass::put_ctree(ostream& fout, int rootidx, int jidx, varclass& var)
// クラスターのケースIDを再帰的に出力する
{
  if (var.cc(rootidx).lidx() < 0) // 2分木の葉(各個体)
  {
    var.cc(rootidx).lno(jidx);
    var.cc(rootidx).rno(jidx);
    fout << "[" << var.cc(rootidx).cid() << "]";
    return jidx + 1;
  }
  jidx = put_ctree(fout, var.cc(rootidx).lidx(), jidx, var);
  jidx = put_ctree(fout, var.cc(rootidx).ridx(), jidx, var);
  return jidx;
}

void clustclass::put_cid(ostream& fout, option& opt, varclass& var)
// クラスターのIDを出力する
{
  for (int i = 0; i < var.nc(); ++i)
  {
    if (var.dc(i, i) > 1.0)     // 2つ以上のケースがあるクラスターだけ出力
    {
      fout << "  融合クラスター[" << var.cc(i).cid() << "]=";
      put_ctree(fout, var.cc(i).ridx(), 0, var);
      fout << "\n";
    }
  }
}

void clustclass::put_dns(option& opt, int rootidx, varclass& var)
// ケースIDとクラスター番号をデータファイル形式で出力する
{
  ofstream fout;
  if (ofopen(fout, (opt.dnsname()).c_str()) == false) return; // ファイルオープン失敗
  idd tid(id.title(), 3);
  tid.assign_item(0, id.item(0));
  ostringstream ostr;
  ostr << "ケース順位(距離:" << opt.disstr() << " 融合手法:" << opt.methstr() << ")";
  tid.assign_item(1, ostr.str());
  ostr.str("");
  ostr << "クラスター番号(距離:" << opt.disstr() << " 融合手法:" << opt.methstr() << ")";
  tid.assign_item(2, ostr.str());
  tid.put(fout, formtag);
  datac cdata(3);
  int gno = 1;
  put_dnsctree(fout, cdata, gno, opt.den(), 0, rootidx, 0, var);
}

int clustclass::put_dnsctree(ofstream& fout, datac& cdata, int& gno, int tlevel0, int tlevel, int rootidx, int jidx, varclass& var)
// ケースIDとクラスター番号を再帰的に出力する
{
  if (var.cc(rootidx).lidx() < 0) // 2分木の葉(各個体)
  {
    ++jidx;
    cdata.assign(0, var.cc(rootidx).cid());
    cdata.assign(1, jidx);
    cdata.assign(2, gno);
    cdata.put(fout);
    return jidx;
  }
  ++tlevel;
  jidx = put_dnsctree(fout, cdata, gno, tlevel0, tlevel, var.cc(rootidx).lidx(), jidx, var);
  if (tlevel <= tlevel0) ++gno;
  jidx = put_dnsctree(fout, cdata, gno, tlevel0, tlevel, var.cc(rootidx).ridx(), jidx, var);
  return jidx;
}

void clustclass::plot_dend(option& opt, int rootidx, varclass& var)
// 樹形図用データを出力する
{
  ofstream fout;
  if (ofopen(fout, (opt.ngpname()).c_str()) == false) return; // ファイルオープン失敗
  fout << "-s1 -dline -w80 -c255,0,0 -zx0," << var.nc() + 1 << ",1\n";
  fout << "%";  puttitle(fout, Title);
  fout << "%\n";
  id.put_title(fout, formngp);
  id.put_stra(fout, formngp);
  fout << "%ケースID:" << id.item(0) << "\n";
  for (size_t i = 1; i < id.ni(); ++i)
  { fout << "%変数x" << setw(3) << i << ":" << id.item(i) << "\n"; }
  fout << "%クラスター間の距離:" << opt.disstr() << "  融合手法:" << opt.methstr() << "\n";
  fout << "%";
  put_ctree(fout, rootidx, 0, var);
  fout << "\n";
  double xc = 0.0, yc = 0.0;
  plot_dtree(fout, xc, yc, rootidx, var);
  fout << xc << " " << yc << "\n";
  fout << xc << " " << yc * 1.1 << endl;
}

void clustclass::plot_dtree(ofstream& fout, double& xc, double& yc, int rootidx, varclass& var) const
// 樹形図の幹データを再帰的に出力する
{
  if (var.cc(rootidx).lidx() < 0) // 2分木の葉(各個体)
  {
    xc = (double)(var.cc(rootidx).rno() + 1);  yc = 0.0;
    return;
  }
  double xl = 0.0, yl = 0.0, xr = 0.0, yr = 0.0;
  plot_dtree(fout, xl, yl, var.cc(rootidx).lidx(), var);
  plot_dtree(fout, xr, yr, var.cc(rootidx).ridx(), var);
  xc = (xl + xr) * 0.5;  yc = var.dc(var.cc(rootidx).lno(), var.cc(rootidx).rno());
  fout << xl << " " << yl << "\n";
  fout << xl << " " << yc << "\n";
  fout << xr << " " << yc << "\n";
  fout << xr << " " << yr << "\n";
  fout << NGPdelm2 << "\n";
}

int clustclass::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) != 0) throw option_error("実行時オプションエラー！");
  return nc;
}

void clustclass::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.nc(), var.mevec(), var.vrvec());
      putmat(cout, "・相関行列", matr, var.nc(), var.smat());
    }
    if (var.resvar() == reserr)
    { cout << "※クラスター分析不可能！\n\n" << endl; }
    else
    {
      cout << "\n・クラスター間の距離:" << opt.disstr() << "\n";
      cout << "  融合手法:" << opt.methstr() << "\n\n";
    }
  }
  if (var.resvar() == reserr) return;
  int rootidx = put_clust(cout, opt, var);
  if (opt.res("res")) cout << endl;
  if (opt.dnsname() != "") put_dns(opt, rootidx, var);
  if (opt.ngpname() != "") plot_dend(opt, rootidx, var);
}

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