C++ ArrayXd::data方法代码示例

void CMATLAB::matGet(string filename,ArrayXd &X,string Xname)
{
    MATFile *pmat=matOpen(filename,string("r"));
    if (!pmat) throw "CMATLAB::importFrom error: matlab.matOpen failed";
    mxArray *pa=matGetVariable(pmat,Xname);
    if (!pa) throw "CMATLAB::importFrom error: matlab.matGetVariable failed";
    int N=mxGetNumberOfElements(pa);
    if (N<=0) throw "CMATLAB::importFrom error: matlab.mxGetNumberOfElements failed";
    X.resize(N);
    memcpy((void *)(mxGetPr(pa)), (void *)X.data(), sizeof((double *)X.data()));
    mxDestroyArray(pa);
    if (!matClose(pmat)) throw "CMATLAB::importFrom error: matlab.matClose failed";
}

    int gesdd(MatrixXd& A, ArrayXd& S, MatrixXd& Vt) {
	int info, mone = -1, m = A.rows(), n = A.cols();
	std::vector<int> iwork(8 * n);
	double wrk;
	if (m < n || S.size() != n || Vt.rows() != n || Vt.cols() != n)
	    throw std::invalid_argument("dimension mismatch in gesvd");
	F77_CALL(dgesdd)("O", &m, &n, A.data(), &m, S.data(), A.data(),
			 &m, Vt.data(), &n, &wrk, &mone, &iwork[0], &info);
	int lwork(wrk);
	std::vector<double> work(lwork);
	F77_CALL(dgesdd)("O", &m, &n, A.data(), &m, S.data(), A.data(),
			 &m, Vt.data(), &n, &work[0], &lwork, &iwork[0], &info);
	return info;
    }

    double glmDist::aic(const ArrayXd& y, const ArrayXd& n, const ArrayXd& mu,
			const ArrayXd& wt, double dev) const {
	int nn = mu.size();
	double ans =
	    ::Rf_asReal(::Rf_eval(::Rf_lang6(as<SEXP>(d_aic),
					     as<SEXP>(NumericVector(y.data(), y.data() + nn)),
					     as<SEXP>(NumericVector(n.data(), n.data() + nn)),
					     as<SEXP>(NumericVector(mu.data(), mu.data() + nn)),
					     as<SEXP>(NumericVector(wt.data(), wt.data() + nn)),
					     PROTECT(::Rf_ScalarReal(dev))), d_rho));
	UNPROTECT(1);
	return ans;
    }

void CMATLAB::matPut(string filename,const ArrayXd &X,string Xname)
{
    MATFile *pmat=matOpen(filename,string("w"));
    if (!pmat) throw "CMATLAB::exportTo error: matOpen failed";
    mxArray *pa=mxCreateDoubleMatrix((int)X.size(),1);
    if (!pa) throw "CMATLAB::exportTo error: mxCreateDoubleMatrix failed";
    memcpy((void *)(mxGetPr(pa)), (void *)X.data(), X.size()*sizeof(double));
    if (!matPutVariable(pmat,string("X"),pa)) throw "CMATLAB::exportTo error: matlab.matPutVariable failed";
    mxDestroyArray(pa);
    if (!matClose(pmat)) throw "CMATLAB::exportTo error: matlab.matClose failed";

}

    void merPredD::updateXwts(const ArrayXd& sqrtXwt) {
        if (d_Xwts.size() != sqrtXwt.size())
            throw invalid_argument("updateXwts: dimension mismatch");
        std::copy(sqrtXwt.data(), sqrtXwt.data() + sqrtXwt.size(), d_Xwts.data());
        if (sqrtXwt.size() == d_V.rows()) { // W is diagonal
            d_V              = d_Xwts.asDiagonal() * d_X;
            for (int j = 0; j < d_N; ++j)
                for (MSpMatrixd::InnerIterator Utj(d_Ut, j), Ztj(d_Zt, j);
                     Utj && Ztj; ++Utj, ++Ztj)
                    Utj.valueRef() = Ztj.value() * d_Xwts.data()[j];
        } else {
            SpMatrixd      W(d_V.rows(), sqrtXwt.size());
            const double *pt = sqrtXwt.data();
            W.reserve(sqrtXwt.size());
            for (Index j = 0; j < W.cols(); ++j, ++pt) {
                W.startVec(j);
                W.insertBack(j % d_V.rows(), j) = *pt;
            }
            W.finalize();
            d_V              = W * d_X;
            SpMatrixd      Ut(d_Zt * W.adjoint());
            if (Ut.cols() != d_Ut.cols())
                throw std::runtime_error("Size mismatch in updateXwts");

            // More complex code to handle the pruning of zeros
            MVec(d_Ut.valuePtr(), d_Ut.nonZeros()).setZero();
            for (int j = 0; j < d_Ut.outerSize(); ++j) {
                MSpMatrixd::InnerIterator lhsIt(d_Ut, j);
                for (SpMatrixd::InnerIterator  rhsIt(Ut, j); rhsIt; ++rhsIt, ++lhsIt) {
                    Index                         k(rhsIt.index());
                    while (lhsIt && lhsIt.index() != k) ++lhsIt;
                    if (lhsIt.index() != k)
                        throw std::runtime_error("Pattern mismatch in updateXwts");
                    lhsIt.valueRef() = rhsIt.value();
                }
            }
        }
        d_VtV.setZero().selfadjointView<Eigen::Upper>().rankUpdate(d_V.adjoint());
        updateL();
    }

    const ArrayXd glmDist::devResid(const ArrayXd &y, const ArrayXd &mu, const ArrayXd &wt) const {
	int n = mu.size();
	return as<ArrayXd>(::Rf_eval(::Rf_lang4(as<SEXP>(d_devRes),
						as<SEXP>(NumericVector(y.data(), y.data() + n)),
						as<SEXP>(NumericVector(mu.data(), mu.data() + n)),
						as<SEXP>(NumericVector(wt.data(), wt.data() + n))
					 ), d_rho));
    }

    const ArrayXd glmDist::variance(const ArrayXd &mu) const {
	return as<ArrayXd>(::Rf_eval(::Rf_lang2(as<SEXP>(d_variance),
						as<SEXP>(Rcpp::NumericVector(mu.data(),
									     mu.data() + mu.size()))
					 ), d_rho));
    }

    const ArrayXd glmLink::muEta(const ArrayXd &eta) const {
	return as<ArrayXd>(::Rf_eval(::Rf_lang2(as<SEXP>(d_muEta),
						as<SEXP>(Rcpp::NumericVector(eta.data(),
									     eta.data() + eta.size()))
					 ), d_rho));
    }

    const ArrayXd glmLink::linkFun(const ArrayXd& mu) const {
	return as<ArrayXd>(::Rf_eval(::Rf_lang2(as<SEXP>(d_linkFun),
						as<SEXP>(Rcpp::NumericVector(mu.data(),
									     mu.data() + mu.size()))
					 ), d_rho));
    }