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

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

MatrixXd NumInt::GuassQaudrature(const int& N, double& a, double& b) {

    int N0=N-1;
    const int N1 = N0+1;
    const int N2 = N0+2;

    VectorXd xu;
    xu.setLinSpaced(N1,-1.0,1.0);


    // Legendre-Gauss-Vandermonde Matrix
    //Matrix<double,N1,N2> L = Matrix<double,N1,N2>::Zero();
    MatrixXd L(N1,N2);
    L = MatrixXd::Zero(N1,N2);

    // Derivative of Legendre-Gauss-Vandermonde Matrix
    //Matrix<double,N1,1> Lp = Matrix<double,N1,1>::Zero();
    VectorXd Lp(N1);
    Lp = VectorXd::Zero(N1);


    VectorXd dum;
    dum.setLinSpaced(N1,0.0,N0);
    ArrayXd y;
    y = cos((2*dum.array()+1)*M_PI/(2*N0+2))+(0.27/N1)*sin(M_PI*xu.array()*N0/N2);

    double deps = std::numeric_limits<double>::epsilon();

    //Initial Guess
    //Array<double,N1,1> y0 = Array<double,N1,1>::Constant(2);
    ArrayXd y0 = ArrayXd::Constant(N1,2);

    while ((y-y0).abs().matrix().maxCoeff() > deps) {


        // L.col(0) = Matrix<double,N1,1>::Constant(1);
        L.col(0) = VectorXd::Constant(N1,1);
        //Lp = Matrix<double,N1,1>::Zero();
        Lp = VectorXd::Zero(N1);

        L.col(1) = y;

        for (int k=1; k!=N1; k++)
        {
            L.col(k+1) = ((2*k+1)*L.col(k).cwiseProduct(y.matrix())-k*L.col(k-1))/(k+1);
        }

        Lp = (N2)*(L.col(N0)-L.col(N1).cwiseProduct(y.matrix())).cwiseQuotient((1-y.square()).matrix());


        y0 = y;
        y = y0-(L.col(N1).cwiseQuotient(Lp)).array();
    }

    // Gauss Points
    //Matrix<double,N1,1> z = ((a*(1-y)+b*(1+y))/2).matrix();
    VectorXd z(N1);
    z = ((a*(1-y)+b*(1+y))/2).matrix();

    // Gauss Weights
    //Matrix<double,N1,1> w;
    VectorXd w(N1);
    w = (b-a)/(((1-y.square()).matrix()).cwiseProduct(Lp.cwiseProduct(Lp))).array()*pow((double)N2/N1,2);

    // Store
    //Matrix<double,N1,2> zw;
    Matrix<double,Dynamic,Dynamic> zw(N1,2);
    zw.col(0)=z;
    zw.col(1)=w;

    return zw;
}

int polyDeg(std::vector<double>& xVec,
			std::vector<double>& yVec,
			std::vector<double>& ysVec)
{
	Map<MatrixXd> x(&xVec[0], xVec.size(), 1);
	Map<MatrixXd> y(&yVec[0], yVec.size(), 1);
	int N = (int)xVec.size();

	double meanY = y.mean();
	MatrixXd ys = MatrixXd::Constant(N,1,meanY);
	// cout << "ys = " << ys << endl;
	// cout << "y = " << y << endl;

	double sum_ys_y_2 = (ys - y).array().square().matrix().sum();
	double tobelog = 2*M_PI*sum_ys_y_2/N;
	double logtobelog = log(tobelog);
	double AIC = 2. + N*(logtobelog +1.) + 4./(N-2.);
	// cout << "AIC= " << AIC << endl;

	MatrixXd p = MatrixXd::Constant(2,2, 0);
	p(0,1) = x.mean();
	// cout << "AIC = " << AIC << endl;
	// cout << "p = " << p << endl;

	MatrixXd PL = MatrixXd::Constant(N,2, 1);
	// cout << "PL = " << PL << endl;

	PL.col(1) = x.array() - p(0,1);
	// cout << "PL = " << PL << endl;

	int n = 1;
	int nit = 0;

	MatrixXd ys1 = ys;
	MatrixXd ys2 = ys;
	MatrixXd ys3 = ys;
	MatrixXd ys4 = ys;

	while (nit<3) {
		int ni = n - 1;
		if (ni>0) {
			ArrayXd x_PL = x.array() * (PL.col(ni).array().square());
			double ppa = x_PL.sum();
			double ppb = PL.col(ni).array().square().sum();
			ArrayXd x_PL2 = x.array() * PL.col(ni-1).array() * PL.col(ni).array();
			p.conservativeResize(p.rows(), p.cols()+1);
			p(0,ni+1) = ppa / ppb;
			p(1,ni+1) = x_PL2.sum() / PL.col(ni-1).array().square().sum();

			PL.conservativeResize(PL.rows(), PL.cols()+1);
			PL.col(ni+1) = (x.array() - p(0,ni+1)) * PL.col(ni).array()
						   - p(1,ni+1)* PL.col(ni-1).array();
			//cout << "PL = " << PL << endl;
		}

		MatrixXd bsxfun_time_y_pl = y.asDiagonal() * PL;
		MatrixXd pl_square = PL.array().square();

		ArrayXd tmp = bsxfun_time_y_pl.colwise().sum().array() / pl_square.colwise().sum().array();
		ys = (PL * tmp.matrix().asDiagonal()).rowwise().sum();

		MatrixXd ys_minus_y = ys - y;
		ArrayXd ys_minus_y_square = ys_minus_y.array().square();
		sum_ys_y_2 = ys_minus_y_square.matrix().sum();

		double log_v = log(2.*M_PI*sum_ys_y_2/N);
		double aic = 2.*(n+1) + N*(log_v + 1.) +
					 2.*(n+1)*(n+2.)/(N-n-2.); // correction for small sample sizes

		// cout << "aic = " << aic << endl;
		if (aic>=AIC) {
			++nit;
		} else {
			nit = 0;
			AIC = aic;
		}
		++n;

		ys1 = ys2;
		ys2 = ys3;
		ys3 = ys4;
		ys4 = ys;

		if (n>=N) {
			ys1 = ys;
			break;
		}
	}

	ysVec.resize(ys1.rows());
	Map<MatrixXd> dest(&ysVec[0],ys1.rows(),1);
	dest = ys1;

	return n - nit - 1;
}

int main() {

    // The eigen approach
    ArrayXd n                = ArrayXd::LinSpaced(N+1,0,N);
    double multiplier        = M_PI/N;
    Array<double, 1, N+1> nT = n.transpose();
    ArrayXd x                = - cos(multiplier*n);
    ArrayXd xsub             = x.middleRows(1, N-1);
    ArrayXd ysub             = (x1-x0)/2*xsub + (x1+x0)/2;

    ArrayXXd T               = cos((acos(x).matrix()*nT.matrix()).array());
    ArrayXXd Tsub            = cos((acos(xsub).matrix()*nT.matrix()).array());
    ArrayXd sqinx            = 1/sqrt(1-xsub*xsub);

    MatrixXd inv1x2          = (sqinx).matrix().asDiagonal();

    // Can't use the following to test elements of inv1x2
    // std::cout << inv1x2(0,0) << "\n";

    MatrixXd Usub            = inv1x2 * sin(((acos(xsub).matrix())*nT.matrix()).array()).matrix();
    MatrixXd dTsub           = Usub*(n.matrix().asDiagonal());
    MatrixXd d2Tsub          = ((sqinx*sqinx).matrix().asDiagonal())*((xsub.matrix().asDiagonal()) * (dTsub.matrix()) - (Tsub.matrix()) * ((n*n).matrix().asDiagonal()));

    MatrixXd d2T(N+1, N+1);
    RowVectorXd a            = (pow((-1),nT))*(nT*nT+1)*(nT*nT)/3;
    RowVectorXd b            = (nT*nT+1)*(nT*nT)/3;
    d2T.middleRows(1,N-1)    = d2Tsub; 
    d2T.row(0)               = a;
    d2T.row(N)               = b;

    MatrixXd D2              = d2T.matrix() * ((T.matrix()).inverse());
    MatrixXd E2              = D2.middleRows(1,N-1).middleCols(1,N-1);
    MatrixXd Y               = ysub.matrix().asDiagonal();
    MatrixXd H               = - (4 / ((x1-x0)*(x1-x0))) * E2 + k*Y;

    Eigen::EigenSolver<Eigen::MatrixXd> HE(H);
    VectorXcd D              = HE.eigenvalues();
    MatrixXcd V              = HE.eigenvectors();
    std::cout << HE.info() << std::endl;

    // Open ofstream
    ofstream Dfile;
    Dfile.open("D-output.txt");

    ofstream Vfile;
    Vfile.open("V-output.txt");

    ofstream V544file;
    V544file.open("V544-output.txt");

    Dfile.precision(15);
    Dfile << D.real() << "\n";

    Vfile.precision(15);
    Vfile << V.real() << "\n";

    V544file.precision(15);

	for(int i = 1; i<N-1; i++)
    {
		V544file << ysub[i-1];
        V544file << " "        << V.col(544).row(i-1).real() << "\n";
	}
    Dfile.close();
    Vfile.close();
	V544file.close();
	system("gnuplot -p plot.gp");
	system("rsvg-convert -w 2000 -o V544-plot.png V544-plot.svg");

}