C++ ArrayXd类代码示例

Vector4d compute_grad(Vector4d beta, VectorXd x, VectorXd y){
  Vector4d grad;
  ArrayXd tmp;
  ArrayXd pred = model_fun(beta, x);

  assert(x.size()==y.size());

  // beta(0)
  tmp = 1 / (1 + exp(-(x.array()-beta(2))/abs(beta(3))));
  tmp *= pred - y.array();
  grad(0) = tmp.sum() / x.size();

  // beta(1)
  tmp = 1 / (1 + exp(-(x.array()-beta(2))/abs(beta(3))));
  tmp = 1 - tmp;
  tmp *= pred - y.array();
  grad(1) = tmp.sum() / x.size();

  // beta(2)
  tmp = -(beta(0)- beta(1)) * (exp((beta(2)-x.array())/abs(beta(3)))/abs(beta(3))) \
        / (1+exp((beta(2)-x.array())/abs(beta(3)))).pow(2);
  tmp *= pred - y.array();
  grad(2) = tmp.sum() / x.size();

  // beta(3)
  tmp = (beta(0) - beta(1)) * (beta(2)-x.array()).pow(2) * sgn(beta(3)) \
        / (abs(beta(3))*pow(beta(3), 2)*(1+exp((beta(2)-x.array())/abs(beta(3)))).pow(2));
  tmp *= pred - y.array();
  grad(3) = tmp.sum() / x.size();

  return grad;
}

    inline ArrayXd lm::Dplus(const ArrayXd& d) {
	ArrayXd   di(d.size());
	double  comp(d.maxCoeff() * threshold());
	for (int j = 0; j < d.size(); ++j) di[j] = (d[j] < comp) ? 0. : 1./d[j];
	m_r          = (di != 0.).count();
	return di;
    }

    //@{
    double                   gammaDist::aic     (const ArrayXd& y, const ArrayXd& n, const ArrayXd& mu,
						 const ArrayXd& wt, double dev) const {
	double   nn(wt.sum());
	double disp(dev/nn);
	double   ans(0), invdisp(1./disp);
	for (int i = 0; i < mu.size(); ++i)
	    ans += wt[i] * ::Rf_dgamma(y[i], invdisp, mu[i] * disp, true);
	return -2. * ans + 2.;
    }

double computeBinWidth(const MatrixXd& positions) {
	// assumes first col of positions corresponds to dominant eigenvect
	ArrayXd firstCol = positions.col(0).array();
	firstCol -= firstCol.mean();
	double SSE = firstCol.matrix().squaredNorm();
	double variance = SSE / firstCol.size();
	double std = sqrt(variance);
	double targetBinsPerStd = (MAX_HASH_VALUE - HASH_VALUE_OFFSET) / TARGET_HASH_SPREAD_STDS;
	return std / targetBinsPerStd;
}

VectorXd probutils::logsumexp (const MatrixXd& X)
{
  const VectorXd mx = X.rowwise().maxCoeff(); // Get max of each row

  // Perform the sum(exp(x - mx)) part
  ArrayXd se = ((X.colwise() - mx).array().exp()).rowwise().sum();

  // return total log(sum(exp(x))) - hoping for return value optimisation
  return (se.log()).matrix() + mx;
}

ArrayXd GoSUM::CModelVariables::hcPoint2ModelPoint(const ArrayXd &x)
{
    if ( x.size()!=mvs.size() )   throw "GoSUM::CModelVariables::hcPoint2ModelPoint error: wrong dimension";
    int j,dim=int(x.size());
    ArrayXd X(dim);

    for ( j=0; j<dim; j++ )
    {   X(j)=mvs[j].generateSampleValue(x(j));  }

    return X;
}

void Functions::modeProfileSinc(RefArrayXd predictions, const RefArrayXd covariates, 
                               const double centroid, const double height, const double resolution)
{
    ArrayXd sincFunctionArgument = Functions::PI*(covariates - centroid)/resolution;
    ArrayXd sincFunction = sincFunctionArgument.sin() / sincFunctionArgument;


    // Multiply the profile by the height in the PSD

    predictions = height*sincFunction.square();
}

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";

}

AFI::AFI(const bool prompt) : SteadyState() {
	if (prompt) cout << "Enter flip-angle (degrees): " << flush;
	double inFlip;
	QI::Read(cin, inFlip);
	m_flip = ArrayXd::Ones(1) * inFlip * M_PI / 180.;
	if (prompt) cout << "Enter TR1 & TR2 (seconds): " << flush;
	ArrayXd temp;
	QI::ReadArray(cin, temp);
	if (temp.rows() != 2)
        QI_EXCEPTION("Must enter 2 TR values.");
	m_TR1 = temp[0]; m_TR2 = temp[1];
}

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;
    }

NOMAD::Point CMADS::ArrayXd2NOMADPoint(const ArrayXd &x)
{
    int i,n=int(x.size());
    NOMAD::Point p(n);
    for ( i=0; i<n; i++ ) p[i]=x(i);
    return p;
}

void GoSUM::CModelVariables::setNTuple(const ArrayXd &X,int _at)
{
    if ( X.size()!=mvs.size() ) throw "GoSUM::CModelVariables::setNTuple error: bad nTupe size";

    int i,N=int(mvs.size());
    for ( i=0; i<N; i++ ) mvs[i].setSampleValue(X(i),_at);
}

    //@{
    double        negativeBinomialDist::aic     (const ArrayXd& y, const ArrayXd& n, const ArrayXd& mu,
						 const ArrayXd& wt, double dev) const {
	return 2. * (wt * (y + d_theta) * (mu + d_theta).log() -
		     y * mu.log() + (y + 1).unaryExpr(Lgamma<double>()) -
		     d_theta * std::log(d_theta) + lgamma(d_theta) -
		     (d_theta + y).unaryExpr(Lgamma<double>())).sum();
    }

void CMT::HistogramNonlinearity::setParameters(const ArrayXd& parameters) {
	if(parameters.size() != mHistogram.size())
		throw Exception("Wrong number of parameters.");

	for(int i = 0; i < mHistogram.size(); ++i)
		mHistogram[i] = parameters[i];
}