text:
//
//  Copyright (c) 2016, Guillaume GODIN
//  All rights reserved.
//
// Redistribution and use in source and binary forms, with or without
// modification, are permitted provided that the following conditions are
// met:
//
//     * Redistributions of source code must retain the above copyright
//       notice, this list of conditions and the following disclaimer.
//     * Redistributions in binary form must reproduce the above
//       copyright notice, this list of conditions and the following
//       disclaimer in the documentation and/or other materials provided
//       with the distribution.
//     * Neither the name of Institue of Cancer Research.
//       nor the names of its contributors may be used to endorse or promote
//       products derived from this software without specific prior written
//       permission.
//
// THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS
// "AS IS" AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT
// LIMITED TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR
// A PARTICULAR PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT
// OWNER OR CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL,
// SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT
// LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE,
// DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY
// THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT
// (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE
// OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
//
// Adding GETAWAY descriptors by Guillaume Godin
// for build & set RDBASE! => export RDBASE=/Users/mbp/Github/rdkit_mine/

#include <GraphMol/RDKitBase.h>
#include <GraphMol/MolTransforms/MolTransforms.h>

#include "GETAWAY.h"
#include "PBF.h"
#include "MolData3Ddescriptors.h"

#include "GraphMol/PartialCharges/GasteigerCharges.h"
#include "GraphMol/PartialCharges/GasteigerParams.h"
#include <Numerics/EigenSolvers/PowerEigenSolver.h>

#include <Numerics/Matrix.h>
#include <Numerics/SquareMatrix.h>
#include <Numerics/SymmMatrix.h>
#include <boost/foreach.hpp>
#include <math.h>
#include <Eigen/Dense>
#include <Eigen/SVD>
#include <iostream>
#include <deque>
#include <Eigen/Core>
#include <Eigen/QR>

using namespace Eigen;

namespace RDKit {

namespace Descriptors {

namespace {

MolData3Ddescriptors moldata3D;

double roundn(double in, int factor) {
  return round(in * pow(10., factor)) / pow(10., factor);
}

double* retreiveMat(MatrixXd matrix) {
  double* arrayd = matrix.data();
  return arrayd;
}

double* retreiveVect(VectorXd matrix) {
  double* arrayd = matrix.data();
  return arrayd;
}

VectorXd getEigenVect(std::vector<double> v) {
  double* varray_ptr = &v[0];
  Map<VectorXd> V(varray_ptr, v.size());
  return V;
}

double round_to_n_digits_(double x, int n) {
  char buff[32];
  sprintf(buff, "%.*g", n, x);  // use g to have significative digits!
  return atof(buff);
}

double round_to_n_digits(double x, int n) {
  double scale = pow(10.0, ceil(log10(abs(x))) + n);
  return round(x * scale) / scale;
}

bool IsClose(double a, double b, int n) {
  bool isclose = false;
  if (fabs(a - b) <= pow(0.1, n) * 1.1) {
    isclose = true;
  }
  return isclose;
}

int countZeros(std::string ta) {
  int nbzero = 0;
  for (unsigned int i = 0; i < ta.length(); i++) {
    if (ta[i] != '0') {
      break;
    }
    nbzero = nbzero + 1;
  }

  return nbzero;
}

bool IsClose2(double a, double b, unsigned int n) {
  bool isclose = false;

  std::string sa, sb;
  std::string ta, tb;
  std::stringstream outa, outb;
  outa << a;
  sa = outa.str();
  outb << b;
  sb = outb.str();
  ta = sa.substr(sa.find(".") + 1);

  tb = sb.substr(sb.find(".") + 1);

  // same number of digits!
  if (countZeros(ta) == countZeros(tb)) {
    if (ta.length() > tb.length()) {
      tb += '0';
    }
    if (ta.length() < tb.length()) {
      ta += '0';
    }
  }
  // std::cout << tb << ":"<<  atoi(tb.c_str()) << ":" << countZeros(ta) <<
  // "\n";
  // std::cout << ta << ":"<<  atoi(ta.c_str()) << ":" << countZeros(ta) <<
  // "\n";

  if (ta.length() == tb.length()) {
    // last digit +/-1 deviation only!)
    if (abs(atoi(ta.c_str()) - atoi(tb.c_str())) < 2) {
      // std::cout << a << "," << b << " ta:"  << ta << "," << "tb:" << tb<<
      // "\n";
      isclose = true;
    }
  }
  return isclose;
}

bool IsClose3(double a, double b, unsigned int n) {
  bool isclose = false;
  std::string sa, sb;
  std::string ta, tb;
  std::stringstream outa, outb;
  outa << a;
  sa = outa.str();
  outb << b;
  sb = outb.str();
  ta = sa.substr(sa.find(".") + 1);
  tb = sb.substr(sb.find(".") + 1);
  // same number of digits!
  if (ta.length() == tb.length()) {
    // last digit +/-1 deviation only!)
    if (abs(atoi(ta.c_str()) - atoi(tb.c_str())) < 2) {
      // std::cout << a << "," << b << " ta:"  << ta << "," << "tb:" << tb<<
      // "\n";
      isclose = true;
    }
  }
  return isclose;
}

// need to clean that code to have always the same output which is not the case
// the classes used the last 2 digit +/- 1 to cluster in the same group (with 3
// digit precision)
// 0.012 and 0.013 are in the same group while 0.014 are not!
std::vector<double> clusterArray(std::vector<double> data, double precision) {
  std::vector<double> Store;

  // sort the input data
  std::sort(data.begin(), data.end());

  // find the difference between each number and its predecessor
  std::vector<double> diffs;

  std::adjacent_difference(data.begin(), data.end(), std::back_inserter(diffs));

  // convert differences into percentage changes
  // std::transform(diffs.begin(), diffs.end(), data.begin(), diffs.begin(),
  //    std::divides<double>());

  int j = 0;
  int count = 0;
  for (unsigned int i = 0; i < data.size(); i++) {
    // std::cout << diffs[i] << ",";
    count++;
    // if a difference exceeds 0.01 <=> 1%, start a new group: if transform is
    // used!
    // use diff not ratio (with 0.003 precision look like it's what it's used in
    // Dragon 6!)
    if (diffs[i] > pow(0.1, precision)) {
      Store.push_back(count);
      count = 0;
      j++;
    }
  }

  return Store;
}

// need to clean that code to have always the same output which is not the case
// the classes used the last 2 digit +/- 1 to cluster in the same group (with 3
// digit precision)
// 0.012 and 0.013 are in the same group while 0.014 are not!
// alternative clustering method not working as in Dragon
std::vector<double> clusterArray2(std::vector<double> data, int precision) {
  std::vector<double> Store;

  // sort the input data descend order!
  std::sort(data.begin(), data.end(), std::greater<double>());
  std::deque<double> B;
  for (unsigned int i = 0; i < data.size(); i++) {
    B.push_back(data[i]);
    // std::cout << data[i] << ",";
  }
  // std::cout << "\n";

  int count = 0;
  while (!B.empty()) {
    double dk = B.front();
    B.pop_front();
    count++;
    // case where the last value is not grouped it goes to the single last group
    if (B.empty()) {
      Store.push_back(count);
    }

    for (unsigned int i = 0; i < B.size(); i++) {
      if (IsClose2(dk, B[i], precision)) {
        count++;
      } else {
        Store.push_back(count);
        for (unsigned int j = 0; j < i; j++) {
          B.pop_front();
        }
        count = 0;
        break;
      }

      if (i == B.size() - 1) {
        Store.push_back(count);
        B.pop_front();
        count = 0;
      }
    }
  }
  return Store;
}

double* GetGeodesicMatrix(double* dist, int lag, int numAtoms) {
  int sizeArray = numAtoms * numAtoms;
  double* Geodesic = new double[sizeArray];
  for (int i = 0; i < sizeArray; i++) {
    if (dist[i] == lag)
      Geodesic[i] = 1;
    else
      Geodesic[i] = 0;
  }

  return Geodesic;
}

JacobiSVD<MatrixXd> getSVD(MatrixXd A) {
  JacobiSVD<MatrixXd> mysvd(A, ComputeThinU | ComputeThinV);
  return mysvd;
}

MatrixXd GetPinv(MatrixXd A) {
  JacobiSVD<MatrixXd> svd = getSVD(A);
  double pinvtoler = 1.e-3;  // choose your tolerance wisely!
  VectorXd vs = svd.singularValues();
  VectorXd vsinv = svd.singularValues();

  for (unsigned int i = 0; i < A.cols(); ++i) {
    if (vs(i) > pinvtoler)
      vsinv(i) = 1.0 / vs(i);
    else
      vsinv(i) = 0.0;
  }

  MatrixXd S = vsinv.asDiagonal();
  MatrixXd Ap = svd.matrixV() * S * svd.matrixU().transpose();
  return Ap;
}

MatrixXd GetCenterMatrix(MatrixXd Mat) {
  VectorXd v = Mat.colwise().mean();
  MatrixXd X = Mat.rowwise() - v.transpose();
  return X;
}

MatrixXd GetHmatrix(MatrixXd X) {
  MatrixXd Weigthed = X.transpose() * X;
  return X * GetPinv(Weigthed) * X.transpose();
}

MatrixXd GetRmatrix(MatrixXd H, MatrixXd DM, int numAtoms) {
  MatrixXd R = MatrixXd::Zero(numAtoms, numAtoms);
  for (int i = 0; i < numAtoms - 1; i++) {
    for (int j = i + 1; j < numAtoms; j++) {
      R(i, j) = sqrt(H(i, i) * H(j, j)) / DM(i, j);
      R(j, i) = R(i, j);
    }
  }

  return R;
}

std::vector<int> GetHeavyList(const ROMol& mol) {
  int numAtoms = mol.getNumAtoms();

  std::vector<int> HeavyList;
  for (int i = 0; i < numAtoms; ++i) {
    const RDKit::Atom* atom = mol.getAtomWithIdx(i);

    if (atom->getAtomicNum() > 1) {
      HeavyList.push_back(1);
    } else
      HeavyList.push_back(0);
  }
  return HeavyList;
}

double* AppendDouble(double* w, double* Append, int length, int pos) {
  for (int i = pos; i < pos + length; i++) {
    w[i] = Append[i - pos];
  }

  return w;
}

double getRCON(MatrixXd R, MatrixXd Adj, int numAtoms) {
  // similar implementation of J. Chem. Inf. Comput. Sci. 2004, 44, 200-209
  // equation 1 or 2 page 201
  // we use instead of atomic absolute values the atomic relative ones as in
  // Dragon
  double RCON = 0.0;
  VectorXd VSR = R.rowwise().sum();
  for (int i = 0; i < numAtoms - 1; ++i) {
    for (int j = i + 1; j < numAtoms; ++j) {
      if (Adj(i, j) > 0) {
        RCON +=
            sqrt(VSR(i) * VSR(j));  // the sqrt is in the sum not over the sum!
      }
    }
  }
  return RCON;
}

double getHATS(double W1, double W2, double H1, double H2) {
  return W1 * H1 * W2 * H2;
}

double getH(double W1, double W2, double H) { return W1 * H * W2; }

double getMax(double* Rk) {
  double RTp = 0;
  for (int j = 0; j < 8; j++) {
    if (Rk[j] > RTp) {
      RTp = Rk[j];
    }
  }
  return RTp;
}

void getGETAWAYDesc(MatrixXd H, MatrixXd R, MatrixXd Adj, int numAtoms,
                    std::vector<int> Heavylist, const ROMol& mol,
                    std::vector<double>& res, int precision) {
  // prepare data for Getaway parameter computation
  // compute parameters

  VectorXd Lev = H.diagonal();
  std::vector<double> heavyLev;

  for (int i = 0; i < numAtoms; i++) {
    if (Heavylist[i] == 1) {
      heavyLev.push_back(round_to_n_digits_(Lev(i), precision));
    }
  }

  std::vector<double> Clus = clusterArray2(heavyLev, precision);

  double numHeavy = heavyLev.size();
  double ITH0 = numHeavy * log(numHeavy) / log(2.0);
  double ITH = ITH0;
  for (unsigned int j = 0; j < Clus.size(); j++) {
    ITH -= Clus[j] * log(Clus[j]) / log(2.0);
  }
  res[0] = roundn(ITH, 3);  // issue somethime with this due to cluster
  double ISH = ITH / ITH0;
  res[1] = roundn(ISH, 3);  // issue somethime with this due to cluster

  // use the PBF to determine 2D vs 3D (with Thresold)
  // determine if this is a plane molecule
  double pbf = RDKit::Descriptors::PBF(mol);
  double D;
  if (pbf < 1.e-5) {
    D = 2.0;
  } else {
    D = 3.0;
  }

  // what about linear molecule ie (D=1) ?
  double HIC = 0.0;
  for (int i = 0; i < numAtoms; i++) {
    HIC -= H(i, i) / D * log(H(i, i) / D) / log(2.);
  }
  res[2] = roundn(HIC, 3);

  double HGM = 1.0;
  for (int i = 0; i < numAtoms; i++) {
    HGM = HGM * H(i, i);
  }
  HGM = 100.0 * pow(HGM, 1.0 / numAtoms);
  res[3] = roundn(HGM, 3);

  double RARS = R.rowwise().sum().sum() / numAtoms;

  JacobiSVD<MatrixXd> mysvd = getSVD(R);

  VectorXd EIG = mysvd.singularValues();

  double rcon = getRCON(R, Adj, numAtoms);

  std::vector<double> wp = moldata3D.GetRelativePol(mol);

  VectorXd Wp = getEigenVect(wp);

  std::vector<double> wm = moldata3D.GetRelativeMW(mol);

  VectorXd Wm = getEigenVect(wm);

  std::vector<double> wi = moldata3D.GetRelativeIonPol(mol);

  VectorXd Wi = getEigenVect(wi);

  std::vector<double> wv = moldata3D.GetRelativeVdW(mol);

  VectorXd Wv = getEigenVect(wv);

  std::vector<double> we = moldata3D.GetRelativeENeg(mol);

  VectorXd We = getEigenVect(we);

  std::vector<double> wu = moldata3D.GetUn(numAtoms);

  VectorXd Wu = getEigenVect(wu);

  std::vector<double> ws = moldata3D.GetIState(mol);

  VectorXd Ws = getEigenVect(ws);

  MatrixXd Bi;
  MatrixXd RBw;
  double HATSu, HATSm, HATSv, HATSe, HATSp, HATSi, HATSs;
  double HATSut, HATSmt, HATSvt, HATSet, HATSpt, HATSit, HATSst;
  double H0ut, H0mt, H0vt, H0et, H0pt, H0it, H0st;
  double R0ut, R0mt, R0vt, R0et, R0pt, R0it, R0st;
  double H0u, H0m, H0v, H0e, H0p, H0i, H0s;
  double R0u, R0m, R0v, R0e, R0p, R0i, R0s;
  double Rkmaxu, Rkmaxm, Rkmaxv, Rkmaxe, Rkmaxp, Rkmaxi, Rkmaxs;
  double tmpu, tmpm, tmpv, tmpe, tmpp, tmpi, tmps;
  double HATSk[7][9];
  double Hk[7][9];
  double Rk[7][8];
  double Rp[7][8];

  double* dist =
      MolOps::getDistanceMat(mol, false);  // need to be be set to false to have
                                           // topological distance not weigthed!

  Map<MatrixXd> D2(dist, numAtoms, numAtoms);

  double Dmax = D2.colwise().maxCoeff().maxCoeff();

  HATSut = 0.0;
  HATSmt = 0.0;
  HATSvt = 0.0;
  HATSet = 0.0;
  HATSpt = 0.0;
  HATSit = 0.0;
  HATSst = 0.0;

  H0ut = 0.0;
  H0mt = 0.0;
  H0vt = 0.0;
  H0et = 0.0;
  H0pt = 0.0;
  H0it = 0.0;
  H0st = 0.0;

  R0ut = 0.0;
  R0mt = 0.0;
  R0vt = 0.0;
  R0et = 0.0;
  R0pt = 0.0;
  R0it = 0.0;
  R0st = 0.0;

  // need to loop other all the D values so <= not <!
  for (int i = 0; i <= Dmax; i++) {
    if (i == 0) {
      Bi = H.diagonal().asDiagonal();
    }

    HATSu = 0.0;
    HATSm = 0.0;
    HATSv = 0.0;
    HATSe = 0.0;
    HATSp = 0.0;
    HATSi = 0.0;
    HATSs = 0.0;

    H0u = 0.0;
    H0m = 0.0;
    H0v = 0.0;
    H0e = 0.0;
    H0p = 0.0;
    H0i = 0.0;
    H0s = 0.0;

    if (i == 0) {  // use Bi
      for (int j = 0; j < numAtoms; j++) {
        for (int k = j; k < numAtoms; k++) {
          if (Bi(j, k) > 0) {
            HATSu += getHATS((double)Wu(j), (double)Wu(j), (double)H(j, j),
                             (double)H(j, j));
            HATSm += getHATS((double)Wm(j), (double)Wm(j), (double)H(j, j),
                             (double)H(j, j));
            HATSv += getHATS((double)Wv(j), (double)Wv(j), (double)H(j, j),
                             (double)H(j, j));
            HATSe += getHATS((double)We(j), (double)We(j), (double)H(j, j),
                             (double)H(j, j));
            HATSp += getHATS((double)Wp(j), (double)Wp(j), (double)H(j, j),
                             (double)H(j, j));
            HATSi += getHATS((double)Wi(j), (double)Wi(j), (double)H(j, j),
                             (double)H(j, j));
            HATSs += getHATS((double)Ws(j), (double)Ws(j), (double)H(j, j),
                             (double)H(j, j));

            if (H(j, k) > 0) {
              H0u += getH((double)Wu(j), (double)Wu(k), (double)H(j, k));
              H0m += getH((double)Wm(j), (double)Wm(k), (double)H(j, k));
              H0v += getH((double)Wv(j), (double)Wv(k), (double)H(j, k));
              H0e += getH((double)We(j), (double)We(k), (double)H(j, k));
              H0p += getH((double)Wp(j), (double)Wp(k), (double)H(j, k));
              H0i += getH((double)Wi(j), (double)Wi(k), (double)H(j, k));
              H0s += getH((double)Ws(j), (double)Ws(k), (double)H(j, k));
            }
          }
        }
      }
    }
    double* Bimat = GetGeodesicMatrix(dist, i, numAtoms);
    Map<MatrixXd> Bj(Bimat, numAtoms, numAtoms);
    if (i > 0 && i < 9) {  // use Bj
      for (int j = 0; j < numAtoms - 1; j++) {
        for (int k = j + 1; k < numAtoms; k++) {
          if (Bj(j, k) == 1) {
            HATSu += getHATS((double)Wu(j), (double)Wu(k), (double)H(j, j),
                             (double)H(k, k));
            HATSm += getHATS((double)Wm(j), (double)Wm(k), (double)H(j, j),
                             (double)H(k, k));
            HATSv += getHATS((double)Wv(j), (double)Wv(k), (double)H(j, j),
                             (double)H(k, k));
            HATSe += getHATS((double)We(j), (double)We(k), (double)H(j, j),
                             (double)H(k, k));
            HATSp += getHATS((double)Wp(j), (double)Wp(k), (double)H(j, j),
                             (double)H(k, k));
            HATSi += getHATS((double)Wi(j), (double)Wi(k), (double)H(j, j),
                             (double)H(k, k));
            HATSs += getHATS((double)Ws(j), (double)Ws(k), (double)H(j, j),
                             (double)H(k, k));

            if (H(j, k) > 0) {
              H0u += getH((double)Wu(j), (double)Wu(k), (double)H(j, k));
              H0m += getH((double)Wm(j), (double)Wm(k), (double)H(j, k));
              H0v += getH((double)Wv(j), (double)Wv(k), (double)H(j, k));
              H0e += getH((double)We(j), (double)We(k), (double)H(j, k));
              H0p += getH((double)Wp(j), (double)Wp(k), (double)H(j, k));
              H0i += getH((double)Wi(j), (double)Wi(k), (double)H(j, k));
              H0s += getH((double)Ws(j), (double)Ws(k), (double)H(j, k));
            }
          }
        }
      }
    }

    if (i >= 9) {  // Totals missing part
      for (int j = 0; j < numAtoms - 1; j++) {
        for (int k = j + 1; k < numAtoms; k++) {
          if (Bj(j, k) == 1) {
            HATSut += 2 * getHATS((double)Wu(j), (double)Wu(k), (double)H(j, j),
                                  (double)H(k, k));
            HATSmt += 2 * getHATS((double)Wm(j), (double)Wm(k), (double)H(j, j),
                                  (double)H(k, k));
            HATSvt += 2 * getHATS((double)Wv(j), (double)Wv(k), (double)H(j, j),
                                  (double)H(k, k));
            HATSet += 2 * getHATS((double)We(j), (double)We(k), (double)H(j, j),
                                  (double)H(k, k));
            HATSpt += 2 * getHATS((double)Wp(j), (double)Wp(k), (double)H(j, j),
                                  (double)H(k, k));
            HATSit += 2 * getHATS((double)Wi(j), (double)Wi(k), (double)H(j, j),
                                  (double)H(k, k));
            HATSst += 2 * getHATS((double)Ws(j), (double)Ws(k), (double)H(j, j),
                                  (double)H(k, k));

            if (H(j, k) > 0) {
              H0ut += 2 * getH((double)Wu(j), (double)Wu(k), (double)H(j, k));
              H0mt += 2 * getH((double)Wm(j), (double)Wm(k), (double)H(j, k));
              H0vt += 2 * getH((double)Wv(j), (double)Wv(k), (double)H(j, k));
              H0et += 2 * getH((double)We(j), (double)We(k), (double)H(j, k));
              H0pt += 2 * getH((double)Wp(j), (double)Wp(k), (double)H(j, k));
              H0it += 2 * getH((double)Wi(j), (double)Wi(k), (double)H(j, k));
              H0st += 2 * getH((double)Ws(j), (double)Ws(k), (double)H(j, k));
            }
          }
        }
      }
    }

    if (i < 9) {
      HATSk[0][i] = HATSu;
      HATSk[1][i] = HATSm;
      HATSk[2][i] = HATSv;
      HATSk[3][i] = HATSe;
      HATSk[4][i] = HATSp;
      HATSk[5][i] = HATSi;
      HATSk[6][i] = HATSs;

      Hk[0][i] = H0u;
      Hk[1][i] = H0m;
      Hk[2][i] = H0v;
      Hk[3][i] = H0e;
      Hk[4][i] = H0p;
      Hk[5][i] = H0i;
      Hk[6][i] = H0s;
    }

    R0u = 0.0;
    R0m = 0.0;
    R0v = 0.0;
    R0e = 0.0;
    R0p = 0.0;
    R0i = 0.0;
    R0s = 0.0;

    Rkmaxu = 0;
    Rkmaxm = 0;
    Rkmaxv = 0;
    Rkmaxe = 0;
    Rkmaxp = 0;
    Rkmaxi = 0;
    Rkmaxs = 0;

    // from 1 to 8
    if (i > 0 && i < 9) {
      for (int j = 0; j < numAtoms - 1; j++) {
        for (int k = j + 1; k < numAtoms; k++) {
          if (Bj(j, k) == 1) {
            tmpu = getH((double)Wu(j), (double)Wu(k),
                        (double)R(j, k));  // Use same function but on all R not
                                           // "H>0" like in the previous loop &
                                           // i>0!
            tmpm = getH((double)Wm(j), (double)Wm(k),
                        (double)R(j, k));  // Use same function but on all R not
                                           // "H>0" like in the previous loop &
                                           // i>0!
            tmpv = getH((double)Wv(j), (double)Wv(k),
                        (double)R(j, k));  // Use same function but on all R not
                                           // "H>0" like in the previous loop &
                                           // i>0!
            tmpe = getH((double)We(j), (double)We(k),
                        (double)R(j, k));  // Use same function but on all R not
                                           // "H>0" like in the previous loop &
                                           // i>0!
            tmpp = getH((double)Wp(j), (double)Wp(k),
                        (double)R(j, k));  // Use same function but on all R not
                                           // "H>0" like in the previous loop &
                                           // i>0!
            tmpi = getH((double)Wi(j), (double)Wi(k),
                        (double)R(j, k));  // Use same function but on all R not
                                           // "H>0" like in the previous loop &
                                           // i>0!
            tmps = getH((double)Ws(j), (double)Ws(k),
                        (double)R(j, k));  // Use same function but on all R not
                                           // "H>0" like in the previous loop &
                                           // i>0!
            R0u += tmpu;
            R0m += tmpm;
            R0v += tmpv;
            R0e += tmpe;
            R0p += tmpp;
            R0i += tmpi;
            R0s += tmps;
            if (tmpu > Rkmaxu) {
              Rkmaxu = tmpu;
            }
            if (tmpm > Rkmaxm) {
              Rkmaxm = tmpm;
            }
            if (tmpv > Rkmaxv) {
              Rkmaxv = tmpv;
            }
            if (tmpe > Rkmaxe) {
              Rkmaxe = tmpe;
            }
            if (tmpp > Rkmaxp) {
              Rkmaxp = tmpp;
            }
            if (tmpi > Rkmaxi) {
              Rkmaxi = tmpi;
            }
            if (tmps > Rkmaxs) {
              Rkmaxs = tmps;
            }
          }
        }
      }

      Rk[0][i - 1] = R0u;
      Rk[1][i - 1] = R0m;
      Rk[2][i - 1] = R0v;
      Rk[3][i - 1] = R0e;
      Rk[4][i - 1] = R0p;
      Rk[5][i - 1] = R0i;
      Rk[6][i - 1] = R0s;

      Rp[0][i - 1] = Rkmaxu;
      Rp[1][i - 1] = Rkmaxm;
      Rp[2][i - 1] = Rkmaxv;
      Rp[3][i - 1] = Rkmaxe;
      Rp[4][i - 1] = Rkmaxp;
      Rp[5][i - 1] = Rkmaxi;
      Rp[6][i - 1] = Rkmaxs;
    }

    if (i >= 9) {
      for (int j = 0; j < numAtoms - 1; j++) {
        for (int k = j + 1; k < numAtoms; k++) {
          if (Bj(j, k) == 1) {
            R0ut += 2 * getH((double)Wu(j), (double)Wu(k),
                             (double)R(j, k));  // Use same function but on all
                                                // R not "H>0" like in the
                                                // previous loop & i>0!
            R0mt += 2 * getH((double)Wm(j), (double)Wm(k),
                             (double)R(j, k));  // Use same function but on all
                                                // R not "H>0" like in the
                                                // previous loop & i>0!
            R0vt += 2 * getH((double)Wv(j), (double)Wv(k),
                             (double)R(j, k));  // Use same function but on all
                                                // R not "H>0" like in the
                                                // previous loop & i>0!
            R0et += 2 * getH((double)We(j), (double)We(k),
                             (double)R(j, k));  // Use same function but on all
                                                // R not "H>0" like in the
                                                // previous loop & i>0!
            R0pt += 2 * getH((double)Wp(j), (double)Wp(k),
                             (double)R(j, k));  // Use same function but on all
                                                // R not "H>0" like in the
                                                // previous loop & i>0!
            R0it += 2 * getH((double)Wi(j), (double)Wi(k),
                             (double)R(j, k));  // Use same function but on all
                                                // R not "H>0" like in the
                                                // previous loop & i>0!
            R0st += 2 * getH((double)Ws(j), (double)Ws(k),
                             (double)R(j, k));  // Use same function but on all
                                                // R not "H>0" like in the
                                                // previous loop & i>0!
          }
        }
      }
    }
    delete[] Bimat;
  }

  // can be column vs row selecgted that can explain the issue!
  double HATSTu = HATSk[0][0] + HATSut;
  double HATSTm = HATSk[1][0] + HATSmt;
  double HATSTv = HATSk[2][0] + HATSvt;
  double HATSTe = HATSk[3][0] + HATSet;
  double HATSTp = HATSk[4][0] + HATSpt;
  double HATSTi = HATSk[5][0] + HATSit;
  double HATSTs = HATSk[6][0] + HATSst;

  for (int ii = 1; ii < 9; ii++) {
    HATSTu += 2 * HATSk[0][ii];
    HATSTm += 2 * HATSk[1][ii];
    HATSTv += 2 * HATSk[2][ii];
    HATSTe += 2 * HATSk[3][ii];
    HATSTp += 2 * HATSk[4][ii];
    HATSTi += 2 * HATSk[5][ii];
    HATSTs += 2 * HATSk[6][ii];
  }

  double HTu = Hk[0][0] + H0ut;
  double HTm = Hk[1][0] + H0mt;
  double HTv = Hk[2][0] + H0vt;
  double HTe = Hk[3][0] + H0et;
  double HTp = Hk[4][0] + H0pt;
  double HTi = Hk[5][0] + H0it;
  double HTs = Hk[6][0] + H0st;

  for (int ii = 1; ii < 9; ii++) {
    HTu += 2.0 * Hk[0][ii];
    HTm += 2.0 * Hk[1][ii];
    HTv += 2.0 * Hk[2][ii];
    HTe += 2.0 * Hk[3][ii];
    HTp += 2.0 * Hk[4][ii];
    HTi += 2.0 * Hk[5][ii];
    HTs += 2.0 * Hk[6][ii];
  }

  // 2*(Rk[1]+Rk[2]+Rk[3]+Rk[4]+Rk[5]+Rk[6]+Rk[7]+Rk[8]) // this is not true
  // this is on all the elemen;

  double RTu = 0.0;
  double RTm = 0.0;
  double RTv = 0.0;
  double RTe = 0.0;
  double RTp = 0.0;
  double RTi = 0.0;
  double RTs = 0.0;

  for (int ii = 0; ii < 8; ii++) {
    RTu += 2.0 * Rk[0][ii];
    RTm += 2.0 * Rk[1][ii];
    RTv += 2.0 * Rk[2][ii];
    RTe += 2.0 * Rk[3][ii];
    RTp += 2.0 * Rk[4][ii];
    RTi += 2.0 * Rk[5][ii];
    RTs += 2.0 * Rk[6][ii];
  }

  double RTMu = getMax(Rp[0]);
  double RTMm = getMax(Rp[1]);
  double RTMv = getMax(Rp[2]);
  double RTMe = getMax(Rp[3]);
  double RTMp = getMax(Rp[4]);
  double RTMi = getMax(Rp[5]);
  double RTMs = getMax(Rp[6]);

  // create the output vector...
  for (int i = 0; i < 9; i++) {
    res[i + 4] = roundn(Hk[0][i], 3);
    res[i + 14] = roundn(HATSk[0][i], 3);
    res[i + 24] = roundn(Hk[1][i], 3);
    res[i + 34] = roundn(HATSk[1][i], 3);
    res[i + 44] = roundn(Hk[2][i], 3);
    res[i + 54] = roundn(HATSk[2][i], 3);
    res[i + 64] = roundn(Hk[3][i], 3);
    res[i + 74] = roundn(HATSk[3][i], 3);
    res[i + 84] = roundn(Hk[4][i], 3);
    res[i + 94] = roundn(HATSk[4][i], 3);
    res[i + 104] = roundn(Hk[5][i], 3);
    res[i + 114] = roundn(HATSk[5][i], 3);
    res[i + 124] = roundn(Hk[6][i], 3);
    res[i + 134] = roundn(HATSk[6][i], 3);
  }

  res[13] = roundn(HTu, 3);
  res[23] = roundn(HATSTu, 3);
  res[33] = roundn(HTm, 3);
  res[43] = roundn(HATSTm, 3);
  res[53] = roundn(HTv, 3);
  res[63] = roundn(HATSTv, 3);
  res[73] = roundn(HTe, 3);
  res[83] = roundn(HATSTe, 3);
  res[93] = roundn(HTp, 3);
  res[103] = roundn(HATSTp, 3);
  res[113] = roundn(HTi, 3);
  res[123] = roundn(HATSTi, 3);
  res[133] = roundn(HTs, 3);
  res[143] = roundn(HATSTs, 3);

  res[144] = roundn(rcon, 3);
  res[145] = roundn(RARS, 3);
  res[146] = roundn(EIG(0), 3);

  for (int i = 0; i < 8; i++) {
    res[i + 147] = roundn(Rk[0][i], 3);
    res[i + 156] = roundn(Rp[0][i], 3);
    res[i + 165] = roundn(Rk[1][i], 3);
    res[i + 174] = roundn(Rp[1][i], 3);
    res[i + 183] = roundn(Rk[2][i], 3);
    res[i + 192] = roundn(Rp[2][i], 3);
    res[i + 201] = roundn(Rk[3][i], 3);
    res[i + 210] = roundn(Rp[3][i], 3);
    res[i + 219] = roundn(Rk[4][i], 3);
    res[i + 228] = roundn(Rp[4][i], 3);
    res[i + 237] = roundn(Rk[5][i], 3);
    res[i + 246] = roundn(Rp[5][i], 3);
    res[i + 255] = roundn(Rk[6][i], 3);
    res[i + 264] = roundn(Rp[6][i], 3);
  }

  res[155] = roundn(RTu + R0ut, 3);
  res[164] = roundn(RTMu, 3);
  res[173] = roundn(RTm + R0mt, 3);
  res[182] = roundn(RTMm, 3);
  res[191] = roundn(RTv + R0vt, 3);
  res[200] = roundn(RTMv, 3);
  res[209] = roundn(RTe + R0et, 3);
  res[218] = roundn(RTMe, 3);
  res[227] = roundn(RTp + R0pt, 3);
  res[236] = roundn(RTMp, 3);
  res[245] = roundn(RTi + R0it, 3);
  res[254] = roundn(RTMi, 3);
  res[263] = roundn(RTs + R0st, 3);
  res[272] = roundn(RTMs, 3);
}

/*
std::vector<std::string>
GETAWAYNAMES={"ITH","ISH","HIC","HGM","H0u","H1u","H2u","H3u","H4u","H5u","H6u","H7u","H8u","HTu",
"HATS0u","HATS1u","HATS2u","HATS3u","HATS4u","HATS5u","HATS6u","HATS7u","HATS8u","HATSu","H0m","H1m","H2m","H3m","H4m","H5m",
"H6m","H7m","H8m","HTm","HATS0m","HATS1m","HATS2m","HATS3m","HATS4m","HATS5m","HATS6m","HATS7m","HATS8m","HATSm","H0v","H1v",
"H2v","H3v","H4v","H5v","H6v","H7v","H8v","HTv","HATS0v","HATS1v","HATS2v","HATS3v","HATS4v","HATS5v","HATS6v","HATS7v","HATS8v",
"HATSv","H0e","H1e","H2e","H3e","H4e","H5e","H6e","H7e","H8e","HTe","HATS0e","HATS1e","HATS2e","HATS3e","HATS4e","HATS5e","HATS6e",
"HATS7e","HATS8e","HATSe","H0p","H1p","H2p","H3p","H4p","H5p","H6p","H7p","H8p","HTp","HATS0p","HATS1p","HATS2p","HATS3p","HATS4p",
"HATS5p","HATS6p","HATS7p","HATS8p","HATSp","H0i","H1i","H2i","H3i","H4i","H5i","H6i","H7i","H8i","HTi","HATS0i","HATS1i","HATS2i",
"HATS3i","HATS4i","HATS5i","HATS6i","HATS7i","HATS8i","HATSi","H0s","H1s","H2s","H3s","H4s","H5s","H6s","H7s","H8s","HTs","HATS0s",
"HATS1s","HATS2s","HATS3s","HATS4s","HATS5s","HATS6s","HATS7s","HATS8s","HATSs","RCON","RARS","REIG","R1u","R2u","R3u","R4u","R5u",
"R6u","R7u","R8u","RTu","R1u+","R2u+","R3u+","R4u+","R5u+","R6u+","R7u+","R8u+","RTu+","R1m","R2m","R3m","R4m","R5m","R6m","R7m",
"R8m","RTm","R1m+","R2m+","R3m+","R4m+","R5m+","R6m+","R7m+","R8m+","RTm+","R1v","R2v","R3v","R4v","R5v","R6v","R7v","R8v","RTv",
"R1v+","R2v+","R3v+","R4v+","R5v+","R6v+","R7v+","R8v+","RTv+","R1e","R2e","R3e","R4e","R5e","R6e","R7e","R8e","RTe","R1e+","R2e+",
"R3e+","R4e+","R5e+","R6e+","R7e+","R8e+","RTe+","R1p","R2p","R3p","R4p","R5p","R6p","R7p","R8p","RTp","R1p+","R2p+","R3p+","R4p+",
"R5p+","R6p+","R7p+","R8p+","RTp+","R1i","R2i","R3i","R4i","R5i","R6i","R7i","R8i","RTi","R1i+","R2i+","R3i+","R4i+","R5i+","R6i+",
"R7i+","R8i+","RTi+","R1s","R2s","R3s","R4s","R5s","R6s","R7s","R8s","RTs","R1s+","R2s+","R3s+","R4s+","R5s+","R6s+","R7s+","R8s+","RTs+"};
 */

void GetGETAWAY(double* dist3D, double* AdjMat, std::vector<double> Vpoints,
                const ROMol& mol, const Conformer& conf,
                std::vector<int> Heavylist, std::vector<double>& res,
                int precision) {
  int numAtoms = conf.getNumAtoms();

  Map<MatrixXd> ADJ(AdjMat, numAtoms, numAtoms);

  Map<MatrixXd> DM(dist3D, numAtoms, numAtoms);

  double* vpoints = &Vpoints[0];

  Map<MatrixXd> matorigin(vpoints, 3, numAtoms);

  MatrixXd MatOrigin = matorigin.transpose();

  MatrixXd Xmean = GetCenterMatrix(MatOrigin);

  MatrixXd H = GetHmatrix(Xmean);

  MatrixXd R = GetRmatrix(H, DM, numAtoms);

  getGETAWAYDesc(H, R, ADJ, numAtoms, Heavylist, mol, res, precision);
}

}  // end of anonymous namespace

void GETAWAY(const ROMol& mol, std::vector<double>& res, int confId,
             int precision) {
  PRECONDITION(mol.getNumConformers() >= 1, "molecule has no conformers")

  int numAtoms = mol.getNumAtoms();

  const Conformer& conf = mol.getConformer(confId);

  std::vector<double> Vpoints(3 * numAtoms);

  for (int i = 0; i < numAtoms; ++i) {
    Vpoints[3 * i] = conf.getAtomPos(i).x;
    Vpoints[3 * i + 1] = conf.getAtomPos(i).y;
    Vpoints[3 * i + 2] = conf.getAtomPos(i).z;
  }

  std::vector<int> Heavylist = GetHeavyList(mol);
  // should be the same as the List size upper!
  // int nHeavyAt= mol.getNumHeavyAtoms();

  double* dist3D = MolOps::get3DDistanceMat(mol, confId);

  double* AdjMat = MolOps::getAdjacencyMatrix(
      mol, false, 0, false, 0);  // false to have only the 1,0 matrix unweighted

  res.clear();
  res.resize(273);

  GetGETAWAY(dist3D, AdjMat, Vpoints, mol, conf, Heavylist, res, precision);
}

}  // end of Descriptors namespace
}  // end of RDKit namespace

max_stars_repo_name: docking-org/rdk

max_issues_repo_name: docking-org/rdk

max_forks_repo_name: docking-org/rdk


text:
#include "Config.hpp"
#include "Util/Unique.hpp"
#include "PopulationGraph.hpp"
#include "Util/Util.hpp"
#include "Util/AssertWithMessage.hpp"
#include "Util/NewExceptionType.hpp"
#include "Util/Strcat.hpp"

#include <unordered_map>
#include <unordered_set>
#include <string>
#include <functional>
#include <algorithm>
#include <array>
#include <utility>
#include <memory>
#include <iterator>
#include <exception>
#include <tuple>
#include <vector>
#include <sstream>

#include <boost/graph/adjacency_list.hpp>
#include <boost/graph/graph_utility.hpp>
#include <boost/graph/filtered_graph.hpp>
#include <boost/graph/connected_components.hpp>
#include <boost/graph/copy.hpp>

namespace pyramid_scheme_simulator {

std::string PopulationGraph::toString() const
{
    std::ostringstream os;

    std::vector<std::pair<Pop, Pop>> allEdges = edges();

    std::multimap<Unique, Unique> edgeMap;
    for(const auto e: allEdges)
    {
        edgeMap.emplace(std::make_pair(e.first->id, e.second->id));
    }

    std::set<Unique> printedSourceVertices;

    //the edge map will sort the uniques
    for(auto e : edgeMap)
    {
        os << e.first << " ---> " << e.second << std::endl;
        printedSourceVertices.emplace(e.first);
    }

    //print the vertices that don't have any edges

    std::set<Unique> allIds;

    {
        const auto allVertices = vertices();
        for(auto v : allVertices)
        {
            allIds.emplace(v->id);
        }
    }

    std::set<Unique> unprintedIds;
    for(auto i : allIds)
    {
        if(printedSourceVertices.find(i) == printedSourceVertices.end())
        {
            unprintedIds.emplace(i);
        }
    }

    for(auto u : unprintedIds)
    {
        os << u << std::endl;
    }

    return os.str();
}

[[noreturn]] void PopulationGraph::throwVertexNotFoundException(const Unique& which)
{
    std::ostringstream os;
    os << "Could not find the index of vertex " << which << std::endl;
    throw PopulationGraph::VertexNotFoundException(os.str());
}

PopulationGraph::BGLPopulationGraph
    PopulationGraph::buildGraph(rd_ptr rd, Config& config)
{

    //check if link chance procs
    auto testEdge = [&rd, &config]() -> bool {
        return Util::sampleFrom(rd,
            config.graphGenerationOptions->linkChance->getOption() / 2.0); };

    const auto graphSize = config.graphGenerationOptions->graphSize.getOption();

    std::unordered_set<UndirectedEdge> edges;
    std::vector<Unique> vertices;

    vertices.reserve(graphSize);

    //generate vertices with random ids
    for(unsigned long i = 0; i < graphSize; i++)
    {
        vertices.push_back(Unique::newUnique());
    }

    //generate edges
    for(auto i : vertices)
    {
        for(auto j : vertices)
        {
            //don't allow edges from a vertex to itself
            if(i != j && testEdge())
            {
                edges.emplace(i, j);
            }
        }
    }

    //build the boost graph
    BGLPopulationGraph g;

    //map a Pop's Unique -> that Pop's vertex descriptor
    std::unordered_map<Unique, BGLPopulationGraph::vertex_descriptor> popDescriptors;

    auto insertPop = [&g, &popDescriptors](Pop p){
        //add it to the graph
        auto v = boost::add_vertex(p, g);
        //then use its id as the key
        popDescriptors.emplace(p->id, v);
    };

    auto startingFunds = config.simulationOptions->startingFunds;

    //add vertices to the boost graph and save the descriptors
    for(auto i : vertices)
    {
        insertPop(std::shared_ptr<CapitalHolder>(
                    config.defaults.mkConsumer(config, rd, i)));
    }

    //add edges using the saved vertex descriptors
    for(auto e : edges)
    {
        //get the vertex descriptors for this edge
        auto u = popDescriptors[e.v1];
        auto v = popDescriptors[e.v2];

        boost::add_edge(u, v, g);
    }

    //TODO: check for disconnected subgraphs

    return g;
}

std::vector<std::unordered_set<PopulationGraph::Pop>>
    PopulationGraph::getDisconnectedSubgraphs(const BGLPopulationGraph& g)
{
    //see http://www.boost.org/doc/libs/1_53_0/libs/graph/example/connected_components.cpp
    //and http://www.boost.org/doc/libs/1_53_0/libs/graph/doc/connected_components.html

    //boost::connected_components assigns every vertex an integer representing its subgraph
    std::vector<int> components;
    components.reserve(boost::num_vertices(g));

    auto numSubgraphs = boost::connected_components(g, &components[0]);

    //initialize the vector of subgraphs and fill it with enough empty sets
    //use a typedef because have to frequently reference ::size_type
    using SubgraphsT = std::vector<std::unordered_set<Pop>>;
    SubgraphsT subgraphs;

    if(numSubgraphs > 0) {
        //if its >0 we can safely cast to an unsigned size type
        subgraphs.reserve(static_cast<SubgraphsT::size_type>(numSubgraphs));

        //impossible to have negative subgraphs
    } else if(numSubgraphs < 0) {
        throw ImplementationException(STRCAT("Error in ", __func__,
                    " call to boost::connected_components, numSubgraphs is < 0"));
    }

    //create empty sets for the subgraphs
    for(int i = 0; i < numSubgraphs; i++)
    {
        //TODO: use fill constructor
        subgraphs.emplace_back();
    }

    //check that every vertex has been assigned to a subgraph
    assert(components.size() == boost::num_vertices(g));

    //assign each vertex to an unordered_set based on the subgraph number it was given
    //by connected_components
    for(std::vector<int>::size_type i = 0;
            i < components.size(); ++i)
    {
        //check each connected component before casting to vector::size_type
        if(components[i] < 0) {
            throw ImplementationException(STRCAT("Error in ", __func__,
                        " components[", i, "] is ", i));
        } else {
            subgraphs[static_cast<SubgraphsT::size_type>(components[i])]
                .insert(g[i]);
        }
    }

    return subgraphs;
}

PopulationGraph::PopulationGraph(Config& config)
{
    graph = buildGraph(config.randomGen, config);
}

PopulationGraph::PopulationGraph(const BGLPopulationGraph& otherGraph)
{
    boost::copy_graph(otherGraph, graph);
}

PopulationGraph::BGLPopulationGraph PopulationGraph::graphFromTuples(
        std::vector<std::pair<Pop, Pop>> tuples)
{
    BGLPopulationGraph g;

    std::unordered_map<Unique, BGLPopulationGraph::vertex_descriptor> popDescriptors;

    //either insert a vertex or return the existing vertex's vertex_descriptor
    const auto insertIfNeeded = [&g, &popDescriptors](Pop& p) {
        const Unique popId = p->id;

        //check if we've already inserted this vertex
        //see http://www.cplusplus.com/reference/unordered_map/unordered_map/count/

        if(popDescriptors.count(popId) == 0) {
            auto vd = boost::add_vertex(p, g);
            //insert the descriptor into the map
            popDescriptors.emplace(popId, vd);

            return vd;
        }
        else {
            return popDescriptors.at(popId);
        }
    };

    for(auto t : tuples)
    {
        auto left = t.first;
        auto right = t.second;

        //get or create the vertex descriptors corresponding to these Pops
        auto leftVd = insertIfNeeded(left);
        auto rightVd = insertIfNeeded(right);

        //use them to create an edge in the graph
        boost::add_edge(leftVd, rightVd, g);
    }

    return g;
}

PopulationGraph::vertices_size_type PopulationGraph::numVertices() const
{
    return boost::num_vertices(graph);
}

PopulationGraph::edges_size_type PopulationGraph::numEdges() const
{
    return boost::num_edges(graph);
}

std::vector<PopulationGraph::Pop> PopulationGraph::vertices() const
{
    //copy all shared_ptr objects into a vector and return it
    std::vector<Pop> allPops;
    allPops.reserve(boost::num_vertices(graph));

    BGLPopulationGraph::vertex_iterator vd, end;
    std::tie(vd, end) = boost::vertices(graph);

    for(; vd != end; ++vd)
    {
        allPops.push_back(graph[*vd]);
    }

    return allPops;
}

std::vector<std::pair<PopulationGraph::Pop, PopulationGraph::Pop>>
    PopulationGraph::edges() const
{
    std::vector<std::pair<PopulationGraph::Pop, PopulationGraph::Pop>>
        allEdges;

    allEdges.reserve(numEdges());

    BGLPopulationGraph::edge_iterator ei, eiEnd;
    std::tie(ei, eiEnd) = boost::edges(graph);

    for(; ei != eiEnd; ++ei)
    {
        allEdges.emplace_back(
            std::make_pair(
                graph[boost::source(*ei, graph)],
                graph[boost::target(*ei, graph)]));
    }

    return allEdges;
}

CapitalHolder& PopulationGraph::findVertexByUnique(const Unique& vert)
{
    return *findVertexPointerByUnique(vert);
}

PopulationGraph::Pop PopulationGraph::findVertexPointerByUnique(const Unique& vert)
{
    BGLPopulationGraph::vertex_iterator begin, end;
    std::tie(begin, end) = boost::vertices(graph);

    //since we're searching for a UUID that has no relation to graph structure
    //there's no reason to use a graph-specific search function
    auto res = std::find_if(begin, end,
        [&vert, this](BGLPopulationGraph::vertex_descriptor vd){
            return this->graph[vd]->id == vert;
        });

    if(res == end)
    {
        throwVertexNotFoundException(vert);
    }
    else
    {
        auto vd = *res;
        Pop foundVert = graph[vd];
        return foundVert;
    }
}

/**
 * can't get this to work with templates
 */
PopulationGraph::vertices_size_type
    PopulationGraph::mutateVerticesOfGraph(MutateVertexFunction mutate,
        BGLPopulationGraph& g)
{
    //counter
    vertices_size_type numMutated = 0;

    BGLPopulationGraph::vertex_iterator vi, viEnd;

    std::tie(vi, viEnd) = boost::vertices(g);

    for(; vi != viEnd; ++vi)
    {
        BGLPopulationGraph::vertex_descriptor vd = *vi;

        //the proper way to access the graph using bundled properties
        //is the [] operator
        //see https://stackoverflow.com/questions/28740974/boost-graph-accessing-properties-through-vertex-descriptor
        g[vd] = mutate(g[vd]);

        numMutated++;
    }

    return numMutated;
}

PopulationGraph::vertices_size_type
    PopulationGraph::mutateVertices(MutateVertexFunction mutate)
{
    return PopulationGraph::mutateVerticesOfGraph(mutate, graph);
}

PopulationGraph::edges_size_type
    PopulationGraph::mutateEdges(MutateEdgeFunction mutate)
{
    return PopulationGraph::mutateEdgesOfGraph(mutate, graph);
}

PopulationGraph::vertices_size_type
    PopulationGraph::mutateVerticesWithPredicate(
        MutateVertexFunction mutate,
        VertexPredicate predicate)
{
    /**
     * orchestrates state to call the VertexPredicate function
     * (i.e. translating between the boost::graph representation
     * of our vertices and CapitalHolder&)
     */
    class VertexPredicateObject
    {
        //TODO: should this be a shared_ptr...?
        //seems like overkill since this is a pointer to a field
        //within the same object this method is running in...
        BGLPopulationGraph* graphPtr = nullptr;

        VertexPredicate pred;

    public:
        //boost filtered_graph docs say the predicate has to have a default constructor
        VertexPredicateObject() = default;

        VertexPredicateObject(
                BGLPopulationGraph& g,
                VertexPredicate _pred
                )
            : graphPtr(&g),
            pred(_pred)
        {}

        /**
         * the docs explicitly say vertex and edge descriptors in the
         * filtered_graph are valid for the original graph object
         */
        bool operator()(const BGLPopulationGraph::vertex_descriptor& vd) const {
            Pop thisVertex = (*graphPtr)[vd];

            //sanity check
            ASSERT_WITH_MESSAGE(thisVertex.get() != nullptr,
                    "The vertices shouldn't have null pointers");

            //dereference the pointer and run the actual predicate
            return pred(*thisVertex);
        }
    } predicateObject(graph, predicate);

    //filter the graph using our predicate
    using FGraph = boost::filtered_graph<BGLPopulationGraph, VertexPredicateObject>;
    FGraph fgraph(graph, predicateObject);

    //***see refactor note at top of mutateVerticesOfGraph***
    vertices_size_type numMutated = 0;

    FGraph::vertex_iterator vi, viEnd;

    std::tie(vi, viEnd) = boost::vertices(fgraph);

    for(; vi != viEnd; ++vi)
    {
        FGraph::vertex_descriptor vd = *vi;

        //the proper way to access the graph using bundled properties
        //is the [] operator
        //see https://stackoverflow.com/questions/28740974/boost-graph-accessing-properties-through-vertex-descriptor
        fgraph[vd] = mutate(fgraph[vd]);
        numMutated++;
    }

    return numMutated;
}

PopulationGraph::edges_size_type
    PopulationGraph::mutateEdgesOfGraph(MutateEdgeFunction mutate,
        BGLPopulationGraph& g)
{
    edges_size_type numMutated = 0;

    BGLPopulationGraph::edge_iterator ei, eiEnd;

    std::tie(ei, eiEnd) = boost::edges(g);

    for(; ei != eiEnd; ++ei)
    {
        BGLPopulationGraph::edge_descriptor ed = *ei;

        Pop s, t;
        std::tie(s, t) = mutate(getVerticesForEdge(ed, g));

        //get the vertex descriptors so we can mutate the vertices in the edge we're
        //currently modifying
        BGLPopulationGraph::vertex_descriptor vdS = boost::source(ed, g);
        BGLPopulationGraph::vertex_descriptor vdT = boost::target(ed, g);

        g[vdS] = s;
        g[vdT] = t;
        numMutated++;
    }

    return numMutated;
}

std::pair<PopulationGraph::Pop, PopulationGraph::Pop>
    PopulationGraph::getVerticesForEdge(
        PopulationGraph::BGLPopulationGraph::edge_descriptor ed,
        BGLPopulationGraph& g)
{
    Pop source = g[boost::source(ed, g)];
    Pop target = g[boost::target(ed, g)];
    return std::make_pair(source, target);
}

PopulationGraph::edges_size_type
    PopulationGraph::mutateEdgesWithPredicate(
        MutateEdgeFunction mutate,
        EdgePredicate predicate)
{
    class EdgePredicateObject
    {
        BGLPopulationGraph* graphPtr;

        EdgePredicate pred;

    public:
        //boost filtered_graph docs say the predicate has to have a default constructor
        EdgePredicateObject() = default;

        EdgePredicateObject(
                BGLPopulationGraph* _graphPtr,
                EdgePredicate _pred
                )
            : graphPtr(_graphPtr),
              pred(_pred)
        {}

        bool operator()(const BGLPopulationGraph::edge_descriptor& ed) const {

            Pop s, t;
            std::tie(s, t) = PopulationGraph::getVerticesForEdge(ed, *graphPtr);

            //sanity check
            ASSERT_WITH_MESSAGE(s.get() != nullptr && t.get() != nullptr,
                    "The edges shouldn't have null pointers");

            //dereference the pointer and run the actual predicate
            return pred(*s, *t);
        }
    } predicateObject(&graph, predicate);

    //filter the graph using our predicate
    using FGraph = boost::filtered_graph<BGLPopulationGraph, EdgePredicateObject>;
    FGraph fgraph(graph, predicateObject);

    edges_size_type numMutated = 0;

    FGraph::edge_iterator ei, eiEnd;

    std::tie(ei, eiEnd) = boost::edges(fgraph);

    for(; ei != eiEnd; ++ei)
    {
        FGraph::edge_descriptor ed = *ei;

        Pop s, t;
        std::tie(s, t) = mutate(getVerticesForEdge(ed, graph));

        //get the vertex descriptors so we can mutate the vertices in the edge we're
        //currently modifying
        FGraph::vertex_descriptor vdS = boost::source(ed, graph);
        FGraph::vertex_descriptor vdT = boost::target(ed, graph);

        fgraph[vdS] = s;
        fgraph[vdT] = t;
        numMutated++;
    }

    return numMutated;
}

void PopulationGraph::auditGraph()
{
    NEW_EXCEPTION_TYPE(AuditGraphException);

    //audit iteration functions
    {
        //setup expected values and checks
        const auto expectedNumVertices = boost::num_vertices(graph);
        const auto expectedNumEdges = boost::num_edges(graph);

        BGLPopulationGraph::vertices_size_type numVertices = 0;
        BGLPopulationGraph::edges_size_type    numEdges = 0;

        /***vertex iteration checks***/
        const auto checkNumVertices =
            [&numVertices, expectedNumVertices](std::string fName) {
                if(numVertices != expectedNumVertices)
                {
                    throw AuditGraphException(STRCAT(
                                "Graph failed audit because expected != actual in ",
                                fName, " expected = ", expectedNumVertices,
                                "actual = ", numVertices));
                }
        };

        const std::function<Pop(Pop)> incNumVertices = [&numVertices](Pop x)
            { numVertices++; return x; };

        //run vertex iteration functions and assert against the results
        forEachVertex(incNumVertices);
        checkNumVertices("forEachVertex");
        numVertices = 0;

        mutateVertices(incNumVertices);
        checkNumVertices("mutateVertices");
        numVertices = 0;

        mutateVerticesWithPredicate(incNumVertices,
                //trivial predicate
                [](const CapitalHolder&) { return true; });
        checkNumVertices("mutateVerticesWithPredicate");
        numVertices=0;
        /*****************************/

        /***edge iteration checks*****/
        const auto checkNumEdges =
            [&numEdges, expectedNumEdges](std::string fName) {
                if(numEdges != expectedNumEdges)
                {
                    throw AuditGraphException(STRCAT(
                                "Graph failed audit because expected != actual in ",
                                fName, " expected = ", expectedNumEdges,
                                "actual = ", numEdges));
                }
            };

        const std::function<std::pair<Pop, Pop>(std::pair<Pop, Pop>)> incNumEdges =
            [&numEdges](std::pair<Pop, Pop> p)
            { numEdges++; return p; };

        //run edge iteration functions and assert against the results
        forEachEdge<std::pair<Pop, Pop>>(incNumEdges);
        checkNumEdges("forEachEdge");
        numEdges = 0;

        mutateEdges(incNumEdges);
        checkNumEdges("mutateEdges");
        numEdges = 0;

        mutateEdgesWithPredicate(incNumEdges,
                //trivial predicate
                [](const CapitalHolder&, const CapitalHolder&) { return true; });
        checkNumEdges("mutateEdgesWithPredicate");
        numEdges=0;

        /*****************************/
    }
}

/***************************************/
/********iterator functions*************/
/***************************************/

PopulationGraph::VertexIt PopulationGraph::vBegin()
{
    return mkTransformationIterator(boost::vertices(graph).first,
            transformVertexIt);
}
PopulationGraph::VertexIt PopulationGraph::vEnd()
{
    return mkTransformationIterator(boost::vertices(graph).second,
            transformVertexIt);
}
PopulationGraph::EdgeIt PopulationGraph::eBegin()
{
    return mkTransformationIterator(boost::edges(graph).first,
            transformEdgeIt);
}
PopulationGraph::EdgeIt PopulationGraph::eEnd()
{
    return mkTransformationIterator(boost::edges(graph).second,
            transformEdgeIt);
}

/***************************************/
/***************************************/
/***************************************/

std::unique_ptr<PopulationGraph>
    PopulationGraph::cloneWithCopyVertexFunction(
            std::function<const Pop(const Pop)> copyVertex)
{
    BGLPopulationGraph newGraph;

    //shallow copy
    boost::copy_graph(graph, newGraph);

    std::unique_ptr<PopulationGraph> newPopulationGraph
        //can't use make_unique because we're calling a protected constructor
        = std::unique_ptr<PopulationGraph>(new PopulationGraph(newGraph));

    //the MutateVertexFunction can change what the boost vertex property
    //(the std::shared_ptr) is pointing to by returning a new shared_ptr object
    newPopulationGraph->mutateVertices(copyVertex);

    return newPopulationGraph;
}

std::unique_ptr<PopulationGraph> PopulationGraph::clone()
{
    const std::function<const Pop(const Pop)> clonePop =
        [](const Pop in) -> const Pop
        {
            return in->clone();
        };
    return cloneWithCopyVertexFunction(clonePop);
}

}

ext: cpp

max_stars_repo_name: tjakway/pyramid-scheme-simulator

max_issues_repo_name: tjakway/pyramid-scheme-simulator

max_forks_repo_name: tjakway/pyramid-scheme-simulator


text:
#include <stan/math/fwd/mat.hpp>
#include <gtest/gtest.h>
#include <boost/random/mersenne_twister.hpp>
#include <boost/math/distributions.hpp>

TEST(ProbDistributionsLkjCorr,fvar_double) {
  using stan::math::fvar;
  boost::random::mt19937 rng;
  int K = 4;
  Eigen::Matrix<fvar<double>,Eigen::Dynamic,Eigen::Dynamic> Sigma(K,K);
  Sigma.setZero();
  Sigma.diagonal().setOnes();
  for (int i = 0; i < K*K; i++)
    Sigma(i).d_ = 1.0;
  fvar<double> eta = stan::math::uniform_rng(0,2,rng);
  fvar<double> f = stan::math::do_lkj_constant(eta, K);
  EXPECT_FLOAT_EQ(f.val_, stan::math::lkj_corr_log(Sigma, eta).val_);
  EXPECT_FLOAT_EQ(2.5177896, stan::math::lkj_corr_log(Sigma, eta).d_);
  eta = 1.0;
  f = stan::math::do_lkj_constant(eta, K);
  EXPECT_FLOAT_EQ(f.val_, stan::math::lkj_corr_log(Sigma, eta).val_);
  EXPECT_FLOAT_EQ(f.d_, stan::math::lkj_corr_log(Sigma, eta).d_);
}

TEST(ProbDistributionsLkjCorrCholesky,fvar_double) {
  using stan::math::fvar;
  boost::random::mt19937 rng;
  int K = 4;
  Eigen::Matrix<fvar<double>,Eigen::Dynamic,Eigen::Dynamic> Sigma(K,K);
  Sigma.setZero();
  Sigma.diagonal().setOnes();
  for (int i = 0; i < K*K; i++)
    Sigma(i).d_ = 1.0;
  fvar<double> eta = stan::math::uniform_rng(0,2,rng);
  fvar<double> f = stan::math::do_lkj_constant(eta, K);
  EXPECT_FLOAT_EQ(f.val_, stan::math::lkj_corr_cholesky_log(Sigma, eta).val_);
  EXPECT_FLOAT_EQ(6.7766843, stan::math::lkj_corr_cholesky_log(Sigma, eta).d_);
  eta = 1.0;
  f = stan::math::do_lkj_constant(eta, K);
  EXPECT_FLOAT_EQ(f.val_, stan::math::lkj_corr_cholesky_log(Sigma, eta).val_);
  EXPECT_FLOAT_EQ(3, stan::math::lkj_corr_cholesky_log(Sigma, eta).d_);
}

TEST(ProbDistributionsLkjCorr,fvar_fvar_double) {
  using stan::math::fvar;
  boost::random::mt19937 rng;
  int K = 4;
  Eigen::Matrix<fvar<fvar<double> >,Eigen::Dynamic,Eigen::Dynamic> Sigma(K,K);
  Sigma.setZero();
  Sigma.diagonal().setOnes();
  for (int i = 0; i < K*K; i++)
    Sigma(i).d_.val_ = 1.0;
  fvar<fvar<double> > eta = stan::math::uniform_rng(0,2,rng);
  fvar<fvar<double> > f = stan::math::do_lkj_constant(eta, K);
  EXPECT_FLOAT_EQ(f.val_.val_, stan::math::lkj_corr_log(Sigma, eta).val_.val_);
  EXPECT_FLOAT_EQ(2.5177896, stan::math::lkj_corr_log(Sigma, eta).d_.val_);
  eta = 1.0;
  f = stan::math::do_lkj_constant(eta, K);
  EXPECT_FLOAT_EQ(f.val_.val_, stan::math::lkj_corr_log(Sigma, eta).val_.val_);
  EXPECT_FLOAT_EQ(f.d_.val_, stan::math::lkj_corr_log(Sigma, eta).d_.val_);
}

TEST(ProbDistributionsLkjCorrCholesky,fvar_fvar_double) {
  using stan::math::fvar;
  boost::random::mt19937 rng;
  int K = 4;
  Eigen::Matrix<fvar<fvar<double> >,Eigen::Dynamic,Eigen::Dynamic> Sigma(K,K);
  Sigma.setZero();
  Sigma.diagonal().setOnes();
  for (int i = 0; i < K*K; i++)
    Sigma(i).d_.val_ = 1.0;
  fvar<fvar<double> > eta = stan::math::uniform_rng(0,2,rng);
  fvar<fvar<double> > f = stan::math::do_lkj_constant(eta, K);
  EXPECT_FLOAT_EQ(f.val_.val_, stan::math::lkj_corr_cholesky_log(Sigma, eta).val_.val_);
  EXPECT_FLOAT_EQ(6.7766843, stan::math::lkj_corr_cholesky_log(Sigma, eta).d_.val_);
  eta = 1.0;
  f = stan::math::do_lkj_constant(eta, K);
  EXPECT_FLOAT_EQ(f.val_.val_, stan::math::lkj_corr_cholesky_log(Sigma, eta).val_.val_);
  EXPECT_FLOAT_EQ(3, stan::math::lkj_corr_cholesky_log(Sigma, eta).d_.val_);
}

max_stars_repo_name: yizhang-cae/torsten

max_issues_repo_name: yizhang-cae/torsten

max_forks_repo_name: yizhang-cae/torsten


text:
#include <boost/lexical_cast.hpp>
#include <string>
#include <iostream>

using namespace std;
using namespace boost;

struct point
{
    int x, y;

    friend bool operator==( point const& l, point const& r )
        { return l.x == r.x && l.y == r.y; }
};

std::ostream& operator<<( std::ostream& out, point const& p )
{
    return out << "(" << p.x << "," << p.y << ")";
}

std::istream& operator>>( std::istream& in, point& result )
{
    char lparen, comma, rparen;
    point p;
    if ( !( in >> lparen >> p.x >> comma >> p.y >> rparen )
            || lparen != '(' || comma != ',' || rparen != ')'
            )
        {
            if ( !( in.exceptions() & std::ios_base::failbit ) )
            {
                in.setstate( std::ios_base::failbit );
                return in;
            }

            throw std::runtime_error( "Error reading Point from stream" );
        }

    result = p;
    return in;
}

int main()
{
    cout << lexical_cast<int>( "42" ) << endl;
    cout << lexical_cast<string>( 42 ) << endl;

    cout << lexical_cast<string>( point{ 10, 21 } ) << endl;
    cout << lexical_cast<point>( "(10,21)" ) << endl;

    point p{ -1, 17 };
    cout << ( lexical_cast<point>( lexical_cast<string>( p ) ) == p ) << endl;
}

max_stars_repo_path: lab28/0 - boost-lexical-cast.cpp

max_issues_repo_path: lab28/0 - boost-lexical-cast.cpp

max_forks_repo_path: lab28/0 - boost-lexical-cast.cpp


text:
#include <iostream>
#include <thread>
#include <pcl/console/parse.h>
#include <pcl/filters/extract_indices.h>
#include <pcl/io/pcd_io.h>
#include <pcl/point_types.h>
#include <pcl/sample_consensus/ransac.h>
#include <pcl/sample_consensus/sac_model_plane.h>
#include <pcl/sample_consensus/sac_model_sphere.h>
#include <pcl/visualization/pcl_visualizer.h>
#include <boost/thread/thread.hpp>

#include "study_sample_consensus.h"

using namespace std::literals::chrono_literals;

static boost::shared_ptr<pcl::visualization::PCLVisualizer>
simpleVis (pcl::PointCloud<pcl::PointXYZ>::ConstPtr cloud)
{
  // --------------------------------------------
  // -----Open 3D viewer and add point cloud-----
  // --------------------------------------------
  boost::shared_ptr<pcl::visualization::PCLVisualizer> viewer (new pcl::visualization::PCLVisualizer ("3D Viewer"));
  viewer->setBackgroundColor (0, 0, 0);
  viewer->addPointCloud<pcl::PointXYZ> (cloud, "sample cloud");
  viewer->setPointCloudRenderingProperties (pcl::visualization::PCL_VISUALIZER_POINT_SIZE, 3, "sample cloud");
  //viewer->addCoordinateSystem (1.0, "global");
  viewer->initCameraParameters ();
  return (viewer);
}

static pcl::visualization::PCLVisualizer::Ptr simpleVis2(pcl::PointCloud<pcl::PointXYZ>::ConstPtr cloud,
                                                         pcl::PointCloud<pcl::PointXYZ>::ConstPtr final = nullptr) {
    // --------------------------------------------
    // -----Open 3D viewer and add point cloud-----
    // --------------------------------------------
    pcl::visualization::PCLVisualizer::Ptr viewer(new pcl::visualization::PCLVisualizer("3D Viewer"));
    viewer->setBackgroundColor(0, 0, 0);
    viewer->addPointCloud<pcl::PointXYZ>(cloud, "sample cloud");
    //像素点大小
    viewer->setPointCloudRenderingProperties(pcl::visualization::PCL_VISUALIZER_POINT_SIZE, 2, "sample cloud");

    if (final != nullptr) {
        pcl::visualization::PointCloudColorHandlerCustom<pcl::PointXYZ> color_handler(final, 255, 0, 0);
        viewer->addPointCloud<pcl::PointXYZ>(final, color_handler, "final cloud");
        viewer->setPointCloudRenderingProperties(pcl::visualization::PCL_VISUALIZER_POINT_SIZE, 4, "final cloud");
    }

    viewer->addCoordinateSystem (1.0, "global");
    viewer->initCameraParameters();
    return (viewer);
}

Study_Sample_Consensus::Study_Sample_Consensus()
{

}

void Study_Sample_Consensus::test_random1(int argc, char ** argv){
    // 初始化点云对象
    pcl::PointCloud<pcl::PointXYZ>::Ptr cloud (new pcl::PointCloud<pcl::PointXYZ>);  //存储源点云
    pcl::PointCloud<pcl::PointXYZ>::Ptr final (new pcl::PointCloud<pcl::PointXYZ>);   //存储提取的局内点

    // 填充点云数据
    cloud->width    = 500;                   //填充点云数目
    cloud->height   = 1;                     //无序点云
    cloud->is_dense = false;                 //是否密集的
    cloud->points.resize (cloud->width * cloud->height);
    for (size_t i = 0; i < cloud->points.size (); ++i)
    {
        if (pcl::console::find_argument (argc, argv, "-s") >= 0 || pcl::console::find_argument (argc, argv, "-sf") >= 0)
        {
            //根据命令行参数用x^2+y^2+Z^2=1设置一部分点云数据，此时点云组成1/4个球体作为内点
            cloud->points[i].x = 1024 * rand () / (RAND_MAX + 1.0);
            cloud->points[i].y = 1024 * rand () / (RAND_MAX + 1.0);
            if (i % 5 == 0) // 500/5 = 100    // 可能会散落在球体外
                cloud->points[i].z = 1024 * rand () / (RAND_MAX + 1.0);   //此处对应的点为局外点
            else if(i % 2 == 0) // (500-100) / 2  ,// 在球体正方向内
                cloud->points[i].z =  sqrt( 1 - (cloud->points[i].x * cloud->points[i].x)
                                            - (cloud->points[i].y * cloud->points[i].y));
            else// 在球体负方向内
                cloud->points[i].z =  - sqrt( 1 - (cloud->points[i].x * cloud->points[i].x)
                                              - (cloud->points[i].y * cloud->points[i].y));
        }
        else
        { //用x+y+z=1设置一部分点云数据，此时地拿云组成的菱形平面作为内点
            cloud->points[i].x = 1024 * rand () / (RAND_MAX + 1.0);
            cloud->points[i].y = 1024 * rand () / (RAND_MAX + 1.0);
            if( i % 2 == 0)
                cloud->points[i].z = 1024 * rand () / (RAND_MAX + 1.0);   //对应的局外点
            else
                cloud->points[i].z = -1 * (cloud->points[i].x + cloud->points[i].y);
        }
    }

    std::vector<int> inliers;  //存储局内点集合的点的索引的向量

    //创建随机采样一致性对象
    pcl::SampleConsensusModelSphere<pcl::PointXYZ>::Ptr
            model_s(new pcl::SampleConsensusModelSphere<pcl::PointXYZ> (cloud));    //针对球模型的对象
    pcl::SampleConsensusModelPlane<pcl::PointXYZ>::Ptr
            model_p (new pcl::SampleConsensusModelPlane<pcl::PointXYZ> (cloud));   //针对平面模型的对象
    if(pcl::console::find_argument (argc, argv, "-f") >= 0)
    {  //根据命令行参数，来随机估算对应平面模型，并存储估计的局内点
        pcl::RandomSampleConsensus<pcl::PointXYZ> ransac (model_p);
        ransac.setDistanceThreshold (.01);           //与平面距离小于0.01 的点称为局内点考虑，注意局内点设置参数
        ransac.computeModel();                       //执行随机参数估计
        ransac.getInliers(inliers);                  //存储估计所得的局内点
    }
    else if (pcl::console::find_argument (argc, argv, "-sf") >= 0 )
    {
        //根据命令行参数  来随机估算对应的圆球模型，存储估计的内点
        pcl::RandomSampleConsensus<pcl::PointXYZ> ransac (model_s);
        ransac.setDistanceThreshold (.01);
        ransac.computeModel();
        ransac.getInliers(inliers);
    }

    // 复制估算模型的所有的局内点到final中
    pcl::copyPointCloud<pcl::PointXYZ>(*cloud, inliers, *final);
    std::cout << *final << std::endl;

    // 创建可视化对象并加入原始点云或者所有的局内点

    boost::shared_ptr<pcl::visualization::PCLVisualizer> viewer;
    if (pcl::console::find_argument (argc, argv, "-f") >= 0 || pcl::console::find_argument (argc, argv, "-sf") >= 0)
        viewer = simpleVis(final);
    else
        viewer = simpleVis(cloud);
    while (!viewer->wasStopped ())
    {
        viewer->spinOnce (100);
        boost::this_thread::sleep (boost::posix_time::microseconds (100000));
    }
   // return 0;
}
void Study_Sample_Consensus::test_random2(int argc, char ** argv){
    // initialize PointClouds
    pcl::PointCloud<pcl::PointXYZ>::Ptr cloud(new pcl::PointCloud<pcl::PointXYZ>);
    pcl::PointCloud<pcl::PointXYZ>::Ptr final(new pcl::PointCloud<pcl::PointXYZ>);

    // populate our PointCloud with points
    cloud->width = 500;
    cloud->height = 1;
    cloud->is_dense = false;
    cloud->points.resize(cloud->width * cloud->height);
    for (std::size_t i = 0; i < cloud->points.size(); ++i) {
        if (pcl::console::find_argument(argc, argv, "-s") >= 0 || pcl::console::find_argument(argc, argv, "-sf") >= 0) {  // -s   -sf
            cloud->points[i].x = 1024 * rand() / (RAND_MAX + 1.0);
            cloud->points[i].y = 1024 * rand() / (RAND_MAX + 1.0);
            if (i % 5 == 0)     // 可能会散落在球体外
                cloud->points[i].z = 1024 * rand() / (RAND_MAX + 1.0);
            else if (i % 2 == 0)// 在球体正方向内
                cloud->points[i].z = sqrt(1 - (cloud->points[i].x * cloud->points[i].x)
                                          - (cloud->points[i].y * cloud->points[i].y));
            else // 在球体负方向内
                cloud->points[i].z = -sqrt(1 - (cloud->points[i].x * cloud->points[i].x)
                                           - (cloud->points[i].y * cloud->points[i].y));
        } else {
            cloud->points[i].x = 1024 * rand() / (RAND_MAX + 1.0);
            cloud->points[i].y = 1024 * rand() / (RAND_MAX + 1.0);
            if (i % 2 == 0)
                cloud->points[i].z = 1024 * rand() / (RAND_MAX + 1.0);
            else
                cloud->points[i].z = -1 * (cloud->points[i].x + cloud->points[i].y);
        }
    }

    std::vector<int> inliers;

    // created RandomSampleConsensus object and compute the appropriated model
    // 球形
    pcl::SampleConsensusModelSphere<pcl::PointXYZ>::Ptr
            model_s(new pcl::SampleConsensusModelSphere<pcl::PointXYZ>(cloud));
    // 平面
    pcl::SampleConsensusModelPlane<pcl::PointXYZ>::Ptr
            model_p(new pcl::SampleConsensusModelPlane<pcl::PointXYZ>(cloud));
    if (pcl::console::find_argument(argc, argv, "-f") >= 0) {  // -f
        pcl::RandomSampleConsensus<pcl::PointXYZ> ransac(model_p);
        ransac.setDistanceThreshold(.01);
        ransac.computeModel();
        ransac.getInliers(inliers);
    } else if (pcl::console::find_argument(argc, argv, "-sf") >= 0) {
        pcl::RandomSampleConsensus<pcl::PointXYZ> ransac(model_s);
        ransac.setDistanceThreshold(.01);
        ransac.computeModel();
        ransac.getInliers(inliers);
    }

    // copies all inliers of the model computed to another PointCloudqq
    // 将cloud中指定索引的点拷贝到final点云中
    pcl::copyPointCloud(*cloud, inliers, *final);

    // creates the visualization object and adds either our original cloud or all of the inliers
    // depending on the command line arguments specified.
    pcl::visualization::PCLVisualizer::Ptr viewer;
    if (pcl::console::find_argument(argc, argv, "-f") >= 0 || pcl::console::find_argument(argc, argv, "-sf") >= 0)
        viewer = simpleVis2(cloud, final);
    else
        viewer = simpleVis2(cloud);
    while (!viewer->wasStopped()) {
        viewer->spinOnce(100);
        std::this_thread::sleep_for(100ms);
    }
}


text:
/**
 *    Copyright (C) 2019-present MongoDB, Inc.
 *
 *    This program is free software: you can redistribute it and/or modify
 *    it under the terms of the Server Side Public License, version 1,
 *    as published by MongoDB, Inc.
 *
 *    This program is distributed in the hope that it will be useful,
 *    but WITHOUT ANY WARRANTY; without even the implied warranty of
 *    MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
 *    Server Side Public License for more details.
 *
 *    You should have received a copy of the Server Side Public License
 *    along with this program. If not, see
 *    <http://www.mongodb.com/licensing/server-side-public-license>.
 *
 *    As a special exception, the copyright holders give permission to link the
 *    code of portions of this program with the OpenSSL library under certain
 *    conditions as described in each individual source file and distribute
 *    linked combinations including the program with the OpenSSL library. You
 *    must comply with the Server Side Public License in all respects for
 *    all of the code used other than as permitted herein. If you modify file(s)
 *    with this exception, you may extend this exception to your version of the
 *    file(s), but you are not obligated to do so. If you do not wish to do so,
 *    delete this exception statement from your version. If you delete this
 *    exception statement from all source files in the program, then also delete
 *    it in the license file.
 */

#include <fstream>
#include <iostream>
#include <memory>

#include <boost/algorithm/string.hpp>
#include <boost/filesystem.hpp>
#include <boost/filesystem/operations.hpp>
#include <boost/format.hpp>
#include <boost/optional/optional_io.hpp>

#include "mongo/bson/json.h"
#include "mongo/client/mongo_uri.h"
#include "mongo/client/sdam/json_test_arg_parser.h"
#include "mongo/client/sdam/sdam_configuration_parameters_gen.h"
#include "mongo/client/sdam/topology_manager.h"
#include "mongo/logv2/log.h"
#include "mongo/stdx/unordered_set.h"
#include "mongo/util/clock_source_mock.h"
#include "mongo/util/options_parser/environment.h"
#include "mongo/util/options_parser/option_section.h"
#include "mongo/util/options_parser/options_parser.h"

#define MONGO_LOGV2_DEFAULT_COMPONENT ::mongo::logv2::LogComponent::kDefault

/**
 * This program runs the Server Discover and Monitoring JSON test files located in
 * the src/monogo/client/sdam/json_tests/sdam_tests sub-directory.
 *
 * The process return code conforms to the UNIX idiom of 0 to indicate success and non-zero to
 * indicate failure. In the case of test failure, the process will return the number of test cases
 * that failed.
 *
 * Example invocation to run all tests:
 *  sdam_json_test --source-dir src/monogo/client/sdam/json_tests/sdam_tests
 *
 * Example invocation to run a single test:
 *  sdam_json_test --source-dir src/monogo/client/sdam/json_tests/sdam_tests --filter
 * normalize_uri_case
 */

namespace fs = boost::filesystem;
namespace moe = mongo::optionenvironment;
using namespace mongo::sdam;

namespace mongo::sdam {

/**
 * This class is responsible for parsing and executing a single 'phase' of the json test
 */
class TestCasePhase {
public:
    TestCasePhase(int phaseNum, MongoURI uri, BSONObj phase) : _testUri(uri), _phaseNum(phaseNum) {
        auto bsonResponses = phase.getField("responses").Array();
        for (auto& response : bsonResponses) {
            const auto pair = response.Array();
            const auto address = HostAndPort(pair[0].String());
            const auto bsonIsMaster = pair[1].Obj();

            if (bsonIsMaster.nFields() == 0) {
                _isMasterResponses.push_back(HelloOutcome(address, BSONObj(), "network error"));
            } else {
                _isMasterResponses.push_back(
                    HelloOutcome(address, bsonIsMaster, duration_cast<HelloRTT>(kLatency)));
            }
        }
        _topologyOutcome = phase["outcome"].Obj();
    }

    // pair of error subject & error description
    using TestPhaseError = std::pair<std::string, std::string>;
    struct PhaseResult {
        std::vector<TestPhaseError> errorDescriptions;
        int phaseNumber;

        bool Success() const {
            return errorDescriptions.size() == 0;
        }
    };
    using PhaseResultPtr = PhaseResult*;

    PhaseResult execute(TopologyManager& topology) const {
        PhaseResult testResult{{}, _phaseNum};

        for (auto response : _isMasterResponses) {
            auto descriptionStr =
                (response.getResponse()) ? response.getResponse()->toString() : "[ Network Error ]";
            LOGV2(20202,
                  "Sending server description",
                  "server"_attr = response.getServer(),
                  "description"_attr = descriptionStr);
            topology.onServerDescription(response);
        }

        LOGV2(20203,
              "TopologyDescription after phase",
              "phaseNumber"_attr = _phaseNum,
              "topologyDescription"_attr = topology.getTopologyDescription()->toString());

        validateServers(
            &testResult, topology.getTopologyDescription(), _topologyOutcome["servers"].Obj());
        validateTopologyDescription(
            &testResult, topology.getTopologyDescription(), _topologyOutcome);

        return testResult;
    }

    int getPhaseNum() const {
        return _phaseNum;
    }

private:
    template <typename T, typename U>
    std::string errorMessageNotEqual(T expected, U actual) const {
        std::stringstream errorMessage;
        errorMessage << "expected '" << actual << "' to equal '" << expected << "'";
        return errorMessage.str();
    }

    std::string serverDescriptionFieldName(const ServerDescriptionPtr serverDescription,
                                           std::string field) const {
        std::stringstream name;
        name << "(" << serverDescription->getAddress() << ") " << field;
        return name.str();
    }

    std::string topologyDescriptionFieldName(std::string field) const {
        std::stringstream name;
        name << "(topologyDescription) " << field;
        return name.str();
    }

    template <typename EVO, typename AV>
    void doValidateServerField(const PhaseResultPtr result,
                               const ServerDescriptionPtr serverDescription,
                               const std::string fieldName,
                               EVO expectedValueObtainer,
                               const AV& actualValue) const {
        const auto expectedValue = expectedValueObtainer();
        if (expectedValue != actualValue) {
            auto errorDescription =
                std::make_pair(serverDescriptionFieldName(serverDescription, fieldName),
                               errorMessageNotEqual(expectedValue, actualValue));
            result->errorDescriptions.push_back(errorDescription);
        }
    }

    void validateServerField(const PhaseResultPtr result,
                             const ServerDescriptionPtr& serverDescription,
                             const BSONElement& expectedField) const {
        const auto serverAddress = serverDescription->getAddress();

        std::string fieldName = expectedField.fieldName();
        if (fieldName == "type") {
            doValidateServerField(result,
                                  serverDescription,
                                  fieldName,
                                  [&]() {
                                      auto status = parseServerType(expectedField.String());
                                      if (!status.isOK()) {
                                          auto errorDescription = std::make_pair(
                                              serverDescriptionFieldName(serverDescription, "type"),
                                              status.getStatus().toString());
                                          result->errorDescriptions.push_back(errorDescription);

                                          // return the actual value since we already have reported
                                          // an error about the parsed server type from the json
                                          // file.
                                          return serverDescription->getType();
                                      }
                                      return status.getValue();
                                  },
                                  serverDescription->getType());

        } else if (fieldName == "setName") {
            doValidateServerField(result,
                                  serverDescription,
                                  fieldName,
                                  [&]() {
                                      boost::optional<std::string> result;
                                      if (expectedField.type() != BSONType::jstNULL) {
                                          result = expectedField.String();
                                      }
                                      return result;
                                  },
                                  serverDescription->getSetName());

        } else if (fieldName == "setVersion") {
            doValidateServerField(result,
                                  serverDescription,
                                  fieldName,
                                  [&]() {
                                      boost::optional<int> result;
                                      if (expectedField.type() != BSONType::jstNULL) {
                                          result = expectedField.numberInt();
                                      }
                                      return result;
                                  },
                                  serverDescription->getElectionIdSetVersionPair().setVersion);

        } else if (fieldName == "electionId") {
            doValidateServerField(result,
                                  serverDescription,
                                  fieldName,
                                  [&]() {
                                      boost::optional<OID> result;
                                      if (expectedField.type() != BSONType::jstNULL) {
                                          result = expectedField.OID();
                                      }
                                      return result;
                                  },
                                  serverDescription->getElectionIdSetVersionPair().electionId);

        } else if (fieldName == "logicalSessionTimeoutMinutes") {
            doValidateServerField(result,
                                  serverDescription,
                                  fieldName,
                                  [&]() {
                                      boost::optional<int> result;
                                      if (expectedField.type() != BSONType::jstNULL) {
                                          result = expectedField.numberInt();
                                      }
                                      return result;
                                  },
                                  serverDescription->getLogicalSessionTimeoutMinutes());

        } else if (fieldName == "minWireVersion") {
            doValidateServerField(result,
                                  serverDescription,
                                  fieldName,
                                  [&]() { return expectedField.numberInt(); },
                                  serverDescription->getMinWireVersion());

        } else if (fieldName == "maxWireVersion") {
            doValidateServerField(result,
                                  serverDescription,
                                  fieldName,
                                  [&]() { return expectedField.numberInt(); },
                                  serverDescription->getMaxWireVersion());

        } else {
            MONGO_UNREACHABLE;
        }
    }

    void validateServers(PhaseResultPtr result,
                         const TopologyDescriptionPtr topologyDescription,
                         const BSONObj bsonServers) const {
        auto actualNumServers = topologyDescription->getServers().size();
        auto expectedNumServers =
            bsonServers.getFieldNames<stdx::unordered_set<std::string>>().size();

        if (actualNumServers != expectedNumServers) {
            std::stringstream errorMessage;
            errorMessage << "expected " << expectedNumServers
                         << " server(s) in topology description. actual was " << actualNumServers
                         << ": ";
            for (const auto& server : topologyDescription->getServers()) {
                errorMessage << server->getAddress() << ", ";
            }
            result->errorDescriptions.push_back(std::make_pair("servers", errorMessage.str()));
        }

        for (const BSONElement& bsonExpectedServer : bsonServers) {
            const auto& serverAddress = HostAndPort(bsonExpectedServer.fieldName());
            const auto& expectedServerDescriptionFields = bsonExpectedServer.Obj();

            const auto& serverDescription = topologyDescription->findServerByAddress(serverAddress);
            if (serverDescription) {
                for (const BSONElement& field : expectedServerDescriptionFields) {
                    validateServerField(result, *serverDescription, field);
                }
            } else {
                std::stringstream errorMessage;
                errorMessage << "could not find server '" << serverAddress
                             << "' in topology description.";
                auto errorDescription = std::make_pair("servers", errorMessage.str());
                result->errorDescriptions.push_back(errorDescription);
            }
        }
    }

    template <typename EVO, typename AV>
    void doValidateTopologyDescriptionField(const PhaseResultPtr result,
                                            const std::string fieldName,
                                            EVO expectedValueObtainer,
                                            const AV& actualValue) const {
        auto expectedValue = expectedValueObtainer();
        if (expectedValue != actualValue) {
            auto errorDescription =
                std::make_pair(topologyDescriptionFieldName(fieldName),
                               errorMessageNotEqual(expectedValue, actualValue));
            result->errorDescriptions.push_back(errorDescription);
        }
    }

    void validateTopologyDescription(PhaseResultPtr result,
                                     const TopologyDescriptionPtr topologyDescription,
                                     const BSONObj bsonTopologyDescription) const {
        {
            constexpr auto fieldName = "topologyType";
            doValidateTopologyDescriptionField(
                result,
                fieldName,
                [&]() { return bsonTopologyDescription[fieldName].String(); },
                toString(topologyDescription->getType()));
        }

        {
            constexpr auto fieldName = "setName";
            doValidateTopologyDescriptionField(result,
                                               fieldName,
                                               [&]() {
                                                   boost::optional<std::string> ret;
                                                   auto bsonField =
                                                       bsonTopologyDescription[fieldName];
                                                   if (!bsonField.isNull()) {
                                                       ret = bsonField.String();
                                                   }
                                                   return ret;
                                               },
                                               topologyDescription->getSetName());
        }

        {
            constexpr auto fieldName = "logicalSessionTimeoutMinutes";
            doValidateTopologyDescriptionField(
                result,
                fieldName,
                [&]() {
                    boost::optional<int> ret;
                    auto bsonField = bsonTopologyDescription[fieldName];
                    if (!bsonField.isNull()) {
                        ret = bsonField.numberInt();
                    }
                    return ret;
                },
                topologyDescription->getLogicalSessionTimeoutMinutes());
        }

        {
            constexpr auto fieldName = "maxSetVersion";
            if (bsonTopologyDescription.hasField(fieldName)) {
                doValidateTopologyDescriptionField(
                    result,
                    fieldName,
                    [&]() {
                        boost::optional<int> ret;
                        auto bsonField = bsonTopologyDescription[fieldName];
                        if (!bsonField.isNull()) {
                            ret = bsonField.numberInt();
                        }
                        return ret;
                    },
                    topologyDescription->getMaxElectionIdSetVersionPair().setVersion);
            }
        }

        {
            constexpr auto fieldName = "maxElectionId";
            if (bsonTopologyDescription.hasField(fieldName)) {
                doValidateTopologyDescriptionField(
                    result,
                    fieldName,
                    [&]() {
                        boost::optional<OID> ret;
                        auto bsonField = bsonTopologyDescription[fieldName];
                        if (!bsonField.isNull()) {
                            ret = bsonField.OID();
                        }
                        return ret;
                    },
                    topologyDescription->getMaxElectionIdSetVersionPair().electionId);
            }
        }

        {
            constexpr auto fieldName = "compatible";
            if (bsonTopologyDescription.hasField(fieldName)) {
                doValidateTopologyDescriptionField(
                    result,
                    fieldName,
                    [&]() { return bsonTopologyDescription[fieldName].Bool(); },
                    topologyDescription->isWireVersionCompatible());
            }
        }
    }

    // the json tests don't actually use this value.
    constexpr static auto kLatency = mongo::Milliseconds(100);

    MongoURI _testUri;
    int _phaseNum;
    std::vector<HelloOutcome> _isMasterResponses;
    BSONObj _topologyOutcome;
};

/**
 * This class is responsible for parsing and executing a single json test file.
 */
class JsonTestCase {
public:
    JsonTestCase(fs::path testFilePath) {
        parseTest(testFilePath);
    }

    struct TestCaseResult {
        std::vector<TestCasePhase::PhaseResult> phaseResults;
        std::string file;
        std::string name;

        bool Success() const {
            return std::all_of(
                phaseResults.begin(),
                phaseResults.end(),
                [](const TestCasePhase::PhaseResult& result) { return result.Success(); });
        }
    };

    TestCaseResult execute() {
        auto config =
            std::make_unique<SdamConfiguration>(getSeedList(),
                                                _initialType,
                                                Milliseconds{kHeartBeatFrequencyMsDefault},
                                                Milliseconds{kConnectTimeoutMsDefault},
                                                Milliseconds{kLocalThresholdMsDefault},
                                                _replicaSetName);

        auto clockSource = std::make_unique<ClockSourceMock>();
        TopologyManagerImpl topology(*config, clockSource.get());

        TestCaseResult result{{}, _testFilePath, _testName};

        for (const auto& testPhase : _testPhases) {
            LOGV2(20204, "### Phase Number ###", "phase"_attr = testPhase.getPhaseNum());
            auto phaseResult = testPhase.execute(topology);
            result.phaseResults.push_back(phaseResult);
            if (!result.Success()) {
                LOGV2(20205, "Phase failed", "phase"_attr = phaseResult.phaseNumber);
                break;
            }
        }

        return result;
    }

    const std::string& Name() const {
        return _testName;
    }

private:
    void parseTest(fs::path testFilePath) {
        _testFilePath = testFilePath.string();
        LOGV2(20207, "### Parsing Test File ###", "testFilePath"_attr = testFilePath.string());
        {
            std::ifstream testFile(_testFilePath);
            std::ostringstream json;
            json << testFile.rdbuf();
            _jsonTest = fromjson(json.str());
        }

        _testName = _jsonTest.getStringField("description").toString();
        _testUri = uassertStatusOK(mongo::MongoURI::parse(_jsonTest["uri"].String()));

        _replicaSetName = _testUri.getOption("replicaSet");
        if (!_replicaSetName) {
            if (_testUri.getServers().size() == 1) {
                _initialType = TopologyType::kSingle;
            } else {
                // We can technically choose either kUnknown or kSharded and be compliant,
                // but it seems that some of the json tests assume kUnknown as the initial state.
                // see: json_tests/sdam_tests/sharded/normalize_uri_case.json
                _initialType = TopologyType::kUnknown;
            }
        } else {
            _initialType = TopologyType::kReplicaSetNoPrimary;
        }

        int phase = 0;
        const std::vector<BSONElement>& bsonPhases = _jsonTest["phases"].Array();
        for (auto bsonPhase : bsonPhases) {
            _testPhases.push_back(TestCasePhase(phase++, _testUri, bsonPhase.Obj()));
        }
    }

    std::vector<HostAndPort> getSeedList() {
        std::vector<HostAndPort> result;
        for (const auto& hostAndPort : _testUri.getServers()) {
            result.push_back(hostAndPort);
        }
        return result;
    }

    BSONObj _jsonTest;
    std::string _testName;
    MongoURI _testUri;
    std::string _testFilePath;
    TopologyType _initialType;
    boost::optional<std::string> _replicaSetName;
    std::vector<TestCasePhase> _testPhases;
};

/**
 * This class runs (potentially) multiple json tests and reports their results.
 */
class SdamJsonTestRunner {
public:
    SdamJsonTestRunner(std::string testDirectory, std::vector<std::string> testFilters)
        : _testFiles(scanTestFiles(testDirectory, testFilters)) {}

    std::vector<JsonTestCase::TestCaseResult> runTests() {
        std::vector<JsonTestCase::TestCaseResult> results;
        const auto testFiles = getTestFiles();
        for (auto jsonTest : testFiles) {
            auto testCase = JsonTestCase(jsonTest);
            try {
                LOGV2(20208, "### Executing Test Case ###", "test"_attr = testCase.Name());
                results.push_back(testCase.execute());
            } catch (const DBException& ex) {
                std::stringstream error;
                error << "Exception while executing " << jsonTest.string() << ": " << ex.toString();
                std::string errorStr = error.str();
                results.push_back(JsonTestCase::TestCaseResult{
                    {TestCasePhase::PhaseResult{{std::make_pair("exception", errorStr)}, 0}},
                    jsonTest.string(),
                    testCase.Name()});
                std::cerr << errorStr;
            }
        }
        return results;
    }

    int report(std::vector<JsonTestCase::TestCaseResult> results) {
        int numTestCases = results.size();
        int numSuccess = 0;
        int numFailed = 0;

        if (std::any_of(
                results.begin(), results.end(), [](const JsonTestCase::TestCaseResult& result) {
                    return !result.Success();
                })) {
            LOGV2(20209, "### Failed Test Results ###");
        }

        for (const auto& result : results) {
            auto file = result.file;
            auto testName = result.name;
            auto phaseResults = result.phaseResults;
            if (result.Success()) {
                ++numSuccess;
            } else {
                LOGV2(20210, "### Test Name ###", "name"_attr = testName);
                LOGV2(20211, "Error in file", "file"_attr = file);
                ++numFailed;
                for (auto phaseResult : phaseResults) {
                    LOGV2(20212, "Phase", "phaseNumber"_attr = phaseResult.phaseNumber);
                    if (!phaseResult.Success()) {
                        for (auto error : phaseResult.errorDescriptions) {
                            LOGV2(20213,
                                  "Errors",
                                  "errorFirst"_attr = error.first,
                                  "errorSecond"_attr = error.second);
                        }
                    }
                }
                LOGV2(20214, "");
            }
        }
        LOGV2(20215,
              "Test cases summary",
              "numTestCases"_attr = numTestCases,
              "numSuccess"_attr = numSuccess,
              "numFailed"_attr = numFailed);

        return numFailed;
    }

    const std::vector<fs::path>& getTestFiles() const {
        return _testFiles;
    }

private:
    std::vector<fs::path> scanTestFiles(std::string testDirectory,
                                        std::vector<std::string> filters) {
        std::vector<fs::path> results;
        for (const auto& entry : fs::recursive_directory_iterator(testDirectory)) {
            if (!fs::is_directory(entry) && matchesFilter(entry, filters)) {
                results.push_back(entry.path());
            }
        }
        return results;
    }

    bool matchesFilter(const fs::directory_entry& entry, std::vector<std::string> filters) {
        const auto filePath = entry.path();
        if (filePath.extension() != ".json") {
            return false;
        }

        if (filters.size() == 0) {
            return true;
        }

        for (const auto& filter : filters) {
            if (filePath.string().find(filter) != std::string::npos) {
                return true;
            } else {
                LOGV2_DEBUG(20216,
                            2,
                            "Test skipped due to filter configuration",
                            "filePath"_attr = filePath.string());
            }
        }

        return false;
    }

    std::vector<fs::path> _testFiles;
};
};  // namespace mongo::sdam

int main(int argc, char* argv[]) {
    ArgParser args(argc, argv);

    ::mongo::logv2::LogManager::global().getGlobalSettings().setMinimumLoggedSeverity(
        MONGO_LOGV2_DEFAULT_COMPONENT, ::mongo::logv2::LogSeverity::Debug(args.Verbose()));

    args.LogParams();

    SdamJsonTestRunner testRunner(args.SourceDirectory(), args.TestFilters());
    return testRunner.report(testRunner.runTests());
}

max_stars_repo_name: benety/mongo

max_issues_repo_name: benety/mongo

max_forks_repo_name: benety/mongo