C++ MatrixXf::ldlt方法代码示例

// Inteprolation de Hardy
void computeHardy(std::vector< std::vector< float >* >* interpoleData, std::vector<std::string>* latitude, std::vector<std::string>* longitude, std::vector<std::string>* data) {
    unsigned int n = data->size();
    Eigen::MatrixXf A = Eigen::MatrixXf(n,n);
    Eigen::VectorXf b = Eigen::VectorXf(data->size());
    Eigen::VectorXf x;

    float lati, longi, latj, longj, latk, longk;
    float minLat, maxLat;
    float minLong, maxLong;
    findExtrema(minLat, maxLat, latitude);
    findExtrema(minLong, maxLong, longitude);

    for (unsigned int i = 0; i < n; ++i) {
        for (unsigned int j = 0; j < n; ++j) {
            lati = atof(latitude->at(i).c_str());
            longi = atof(longitude->at(i).c_str());
            latj = atof(latitude->at(j).c_str());
            longj = atof(longitude->at(j).c_str());
            A(i,j) = sqrt(R + pow(distance(longi,lati,longj,latj),2));
        }
    }

    for (unsigned int i = 0; i < data->size(); ++i) {
        b(i) = atof(data->at(i).c_str());
    }

    x = A.ldlt().solve(b);

    for (int i = 0; i < resolution; ++i) {
        for (int j = 0; j < resolution; ++j) {
            lati = minLat + ((double(i)/(resolution-1)) * (maxLat - minLat));
            longi = minLong + ((double(j)/(resolution-1)) * (maxLong - minLong));
            float eval = 0.0;
            for (unsigned int k = 0; k < n; ++k) {
                latk = atof(latitude->at(k).c_str());
                longk = atof(longitude->at(k).c_str());
                eval += x(k) * sqrt(R + pow(distance(lati, longi, latk, longk),2));
            }
            interpoleData->at(i)->at(j) = eval;
        }
    }
}

bool LogisticRegressionVT::LogisticVTImpl::TestCovariate(const Matrix& Xnull,
                                                         const Vector& yVec,
                                                         const Matrix& Xcol) {
  // const int n = Xnull.rows;
  const int d = Xnull.cols;
  const int k = Xcol.cols;
  copy(Xnull, &cov);  // Z = n by d = [z_1^T; z_2^T; ...]
  copy(Xcol, &geno);  // S = n by k = [S_1^T, S_2^T, ...]
  copy(yVec, &y);
  copy(this->null.GetPredicted(), &res);  // n by 1
  v = res.array() * (1. - res.array());
  res = y - res;
  Eigen::MatrixXf vsz(d, k);  // \sum_i v_i S_ki Z_i = n by d matrix
  Eigen::MatrixXf tmp;

  // calculate U and V
  const Eigen::MatrixXf& S = geno;
  const Eigen::MatrixXf& Z = cov;

  // U = (S.transpose() * (res.asDiagonal())).rowsum();
  U = res.transpose() * S;  // 1 by k matrix
                            // for (int i = 0; i < d; ++i) {
  //   vsz.col(i) = (Z * (v.array() *
  //   S.col().array()).matrix().asDiagonal()).rowsum();
  // }
  vsz = cov.transpose() * v.asDiagonal() * S;  // vsz: d by k matrix
  // const double zz = (v.array() *
  // (Z.array().square().rowise().sum()).array()).sum();
  Eigen::MatrixXf zz =
      cov.transpose() * v.asDiagonal() * cov;  // zz: d by d matrix

  // V.size(k, 1);
  //   V(i, 1) = (v.array() * (S.col(i).array().square())).sum() -
  //   vsz.row(i).transpose() * vsz.row(i) /   zz;
  // }
  V = (v.asDiagonal() * (S.array().square().matrix()))
          .colwise()
          .sum();  // - // 1 by k
  tmp = ((vsz).array() * (zz.ldlt().solve(vsz)).array()).colwise().sum();
  V -= tmp;

  // V = (v.asDiagonal() * (S.array().square().matrix())).colwise().sum() -
  // ((vsz).array() * (zz.ldlt().solve(vsz)).array()).colwise().sum();

  // Uk is n by k matrix
  // Uk.size(n, k);
  // for (int i = 0; i < k; ++i) {
  //   Uk.col(i) = res * (S.col(i) - vsz.col(i).transpose()) * Z.col(i) / zz;
  // }
  Uk = res.asDiagonal() * (S - Z * zz.ldlt().solve(vsz));
  // Vkk.size(k, k);
  // for (int i = 0; i < k; ++i) {
  //   for (int j = 0; j <= 1; ++j) {
  //     Vkk(i, j) = Uk.col(i) .transpose() * Uk.col(j);
  //   }
  //   if (i != j) {
  //     Vkk(j, i) = Vkk(i, j);
  //   } else {
  //     if (Vkk(i,i) == 0.0) {
  //       return false; // variance term should be larger than zero.
  //     }
  //   }
  //  }
  Vkk = Uk.transpose() * Uk;

  Eigen::MatrixXf t = U.array() / V.array().sqrt();
  Eigen::RowVectorXf tmp2 = t.row(0).cwiseAbs();
  tmp2.maxCoeff(&maxIndex);

  rep(-tmp2(maxIndex), k, &lower);
  rep(tmp2(maxIndex), k, &upper);
  makeCov(Vkk);
  if (mvnorm.getBandProbFromCor(k, (double*)lower.data(), (double*)upper.data(),
                                (double*)cor.data(), &pvalue)) {
    fprintf(stderr, "Cannot get MVN pvalue.\n");
    return false;
  }

  copy(U, &retU);
  copy(V, &retV);
  copy(t, &retT);
  copy(Vkk, &retCov);
  pvalue = 1.0 - pvalue;
  return true;
};