C++ VectorXd::cwiseAbs方法代码示例

double EEMS::test_prior(const MatrixXd &mSeeds, const VectorXd &mEffcts, const double mrateMu,
			const MatrixXd &qSeeds, const VectorXd &qEffcts,
			const double df, const double sigma2, const double mrateS2, const double qrateS2) const {
  bool inrange = true;
  int qtiles = qEffcts.size();
  int mtiles = mEffcts.size();
  // First check that all parameters fall into their support range
  for ( int i = 0 ; i < qtiles ; i++ ) {
    if (!habitat.in_point(qSeeds(i,0),qSeeds(i,1))) { inrange = false; }
  }
  for ( int i = 0 ; i < mtiles ; i++ ) {
    if (!habitat.in_point(mSeeds(i,0),mSeeds(i,1))) { inrange = false; }
  }
  if (qEffcts.cwiseAbs().minCoeff()>params.qEffctHalfInterval) { inrange = false; }
  if (mEffcts.cwiseAbs().minCoeff()>params.mEffctHalfInterval) { inrange = false; }
  if (abs(mrateMu)>params.mrateMuHalfInterval) { inrange = false; }
  if ((df<params.dfmin) || (df>params.dfmax)) { inrange = false; }
  if (!inrange) { return (-Inf); }
  // Then compute the prior, on the log scale
  double logpi = - log(df)
    + dnegbinln(mtiles,params.negBiSize,params.negBiProb)
    + dnegbinln(qtiles,params.negBiSize,params.negBiProb)
    + dinvgamln(mrateS2,params.mrateShape_2,params.mrateScale_2)
    + dinvgamln(qrateS2,params.qrateShape_2,params.qrateScale_2)
    + dmvnormln(mEffcts,VectorXd::Zero(mtiles),mrateS2*MatrixXd::Identity(mtiles,mtiles))
    + dmvnormln(qEffcts,VectorXd::Zero(qtiles),qrateS2*MatrixXd::Identity(qtiles,qtiles))
    + dinvgamln(sigma2,params.sigmaShape_2,params.sigmaScale_2);
  return (logpi);
}

double computeMerit(double& delta, StdVectorX& X, StdVectorU& U, double penalty_coeff) {
	double merit = computeObjective(delta, X, U);
	for(int t = 0; t < T-1; ++t) {
		VectorXd hval = dynamics_difference(continuous_cartpole_dynamics, X[t], X[t+1], U[t], delta);
		merit += penalty_coeff*(hval.cwiseAbs()).sum();
	}
	return merit;
}

void CostCalculator_eigen::calculateCost(const real_1d_array &xWeight, double &func, real_1d_array &grad) {
    for (int i = 0; i < variableNumber_; ++i) xReal_(i) = xWeight[i];

    // risk grad
    VectorXd riskGrad = varMatrix_ * xReal_;

    // weight change
    VectorXd weightChange = xReal_ - currentWeight_;

    func = 0.5 * xReal_.dot(riskGrad) + weightChange.cwiseAbs().dot(tradingCost_) - expectReturn_.dot(xReal_);
    VectorXd tmp = riskGrad + weightChange.cwiseSign().cwiseProduct(tradingCost_) - expectReturn_;

    for (int i = 0; i < variableNumber_; ++i) grad(i) = tmp(i);
}

//.........这里部分代码省略.........


    // Check the results of the Cox fit; mirrored from the same checks
    // in coxph.fit.S and coxph.R from the R survival package.

    // A vector to indicate for which covariates the betas/se_betas
    // should be set to NAN.
    std::vector<bool> setToNAN = std::vector<bool>(X.nrow, false);

    // Based on coxph.fit.S lines with 'which.sing' and coxph.R line
    // with if(any(is.NA(coefficients))). These lines set coefficients
    // to NA if flag < nvar (with nvar = ncol(x)) and MAXITER >
    // 0. coxph.R then checks for any NAs in the coefficients and
    // outputs the warning message if NAs were found.
    if (flag < X.nrow)
    {
        int which_sing = 0;
        MatrixXd imateigen = imat.data;
        VectorXd imatdiag = imateigen.diagonal();

        for (int i = 0; i < imatdiag.size(); i++)
        {
            if (imatdiag[i] == 0)
            {
                which_sing = i;
                setToNAN[which_sing] = true;
                if (i != 0) {
                    // Don't warn for i=0 to exclude the beta
                    // coefficient for the (constant) mean from the
                    // check. For Cox regression the constant terms
                    // are ignored. However, we leave it in the
                    // calculations because otherwise the null model
                    // calculation will fail in case there are no
                    // other covariates than the SNP.
                    std::cerr << "Warning for " << snpinfo.name[cursnp]
                              << ", model " << modelNames[model]
                              << ": X matrix deemed to be singular (variable "
                              << which_sing + 1 << ")" << std::endl;
                }
            }
        }
    }

    if (niter >= MAXITER)
    {
        cerr << "Warning for " << snpinfo.name[cursnp]
             << ", model " << modelNames[model]
             << ": nr of iterations > the maximum (" << MAXITER << "): "
             << niter << endl;
    }

    if (flag == 1000)
    {
        cerr << "Warning for " << snpinfo.name[cursnp]
             << ", model " << modelNames[model]
             << ": Cox regression ran out of iterations and did not converge,"
             << " setting beta and se to 'NaN'\n";

        std::fill(setToNAN.begin(), setToNAN.end(), true);
    } else {
        VectorXd ueigen = u.data;
        MatrixXd imateigen = imat.data;
        VectorXd infs = ueigen.transpose() * imateigen;
        infs = infs.cwiseAbs();
        VectorXd betaeigen = beta.data;

        assert(betaeigen.size() == infs.size());

        // We check the beta's for all coefficients
        // (incl. covariates), maybe stick to only checking the SNP
        // coefficient?
        for (int i = 0; i < infs.size(); i++) {
            if (infs[i] > EPS &&
                infs[i] > sqrt(EPS) * abs(betaeigen[i])) {
                setToNAN[i] = true;
                cerr << "Warning for " << snpinfo.name[cursnp]
                     << ", model " << modelNames[model]
                     << ": beta for covariate " << i + 1 << " may be infinite,"
                     << " setting beta and se to 'NaN'\n";
            }
        }
    }

    for (int i = 0; i < X.nrow; i++)
    {
        if (setToNAN[i])
        {
            // Cox regression failed somewhere, set results to NAN for
            // this X row (covariate or SNP)
            sebeta[i] = NAN;
            beta[i]   = NAN;
            loglik    = NAN;
        } else {
            sebeta[i] = sqrt(imat.get(i, i));
            loglik = loglik_int[1];
        }
    }

    delete[] work;
}

//.........这里部分代码省略.........

                    // Normalize new message
                    if (newMessage.sum())
                        newMessage = newMessage / newMessage.sum();

                    //cout << "New message: " << endl << newMessage << endl;
                }

                //
                // Set the message!
                //

                VectorXd smoothedOldMessage(newMessage.rows());
                smoothedOldMessage.setZero();

                double smoothing = options.particularD["smoothing"];

                if ( smoothing != 0 )
                    if ( nodeID == ID1 )
                        newMessage = newMessage + (1-smoothing) * messages[ edgeIndex ][0];
                    else
                        newMessage = newMessage + (1-smoothing) * messages[ edgeIndex ][1];

                //cout << "New message:" << endl << newMessage << endl << "Smoothed" << endl << smoothedOldMessage << endl;

                // If residual belief propagation is activated, just check if the
                // newMessage is the one with the higest residual till the
                // moment. Otherwise, set the new message as the current one
                if ( options.particularS["order"] == "RBP" )
                {                    
                    if ( nodeID == ID1 )
                    {
                        VectorXd differences = messages[edgeIndex][0] - newMessage;
                        double difference = differences.cwiseAbs().sum();

                        if ( difference > maxDifference )
                        {
                            from1to2 = true;
                            edgeWithMaxDiffIndex = edgeIndex;
                            maxDifference = difference;
                            associatedMessage = newMessage;
                        }
                    }
                    else
                    {
                        VectorXd differences = messages[edgeIndex][1] - newMessage;
                        double difference = differences.cwiseAbs().sum();

                        if ( difference > maxDifference )
                        {
                            from1to2 = false;
                            edgeWithMaxDiffIndex = edgeIndex;
                            maxDifference = difference;
                            associatedMessage = newMessage;
                        }
                    }
                }
                else
                {
//                    cout << newMessage.cols() << " " << newMessage.rows() << endl;
//                    cout << "edgeIndex" << edgeIndex << endl;
                    if ( nodeID == ID1 )
                    {
//                        cout << messages[ edgeIndex ][0].cols() << " " << messages[ edgeIndex ][0].rows() << endl;
                        messages[ edgeIndex ][0] = newMessage;
                    }