C++ MatrixView::subMatrix方法代码示例

    void SBBox::SBBoxImpl::fillXValue(tmv::MatrixView<double> val,
                                      double x0, double dx, int izero,
                                      double y0, double dy, int jzero) const
    {
        dbg<<"SBBox fillXValue\n";
        dbg<<"x = "<<x0<<" + i * "<<dx<<", izero = "<<izero<<std::endl;
        dbg<<"y = "<<y0<<" + j * "<<dy<<", jzero = "<<jzero<<std::endl;

        assert(val.stepi() == 1);
        const int m = val.colsize();
        const int n = val.rowsize();
        typedef tmv::VIt<double,1,tmv::NonConj> It;

        // It will be useful to do everything in units of dx,dy
        x0 /= dx;
        double wo2 = _wo2 / std::abs(dx);
        y0 /= dy;
        double ho2 = _ho2 / std::abs(dy);
        xdbg<<"x0,y0 -> "<<x0<<','<<y0<<std::endl;
        xdbg<<"width,height -> "<<wo2*2.<<','<<ho2*2.<<std::endl;

        // Start by setting everything to zero
        val.setZero();

        // Then fill the interior with _norm:
        // Fill pixels where:
        //     x0 + ix >= -width/2
        //     x0 + ix < width/2
        //     y0 + iy >= -width/2
        //     y0 + iy < width/2

        int ix1 = std::max(0, int(std::ceil(-wo2 - x0)));
        int ix2 = std::min(m, int(std::ceil(wo2 - x0)));
        int iy1 = std::max(0, int(std::ceil(-ho2 - y0)));
        int iy2 = std::min(n, int(std::ceil(ho2 - y0)));

        if (ix1 < ix2 && iy1 < iy2)
            val.subMatrix(ix1,ix2,iy1,iy2).setAllTo(_norm);

#if 0
        // We used to implement this by making the pixels that cross the edge have a
        // fractional flux value appropriate for the fraction of the box that goes through
        // each pixel.  However, this isn't actually correct.  SBProfile objects are always
        // rendered as the local surface brightness at the center of the pixel.  To get
        // the right flux, you need to convolve by a Pixel.  So if someone renders a Box
        // without convolving by a pixel, it is inconsistent to do this differently than we
        // do all the other SBProfile types.  However, since it was an involved calculation
        // and someone might actually care to resurrect it in a different guise at some point,
        // I'm leaving it here, just commented out.

        // We need to make sure the pixels where the edges of the box fall only get
        // a fraction of the flux.
        //
        // We divide up the range into 3 sections in x:
        //    left of the box where val = 0
        //    in the box where val = _norm
        //    right of the box where val = 0 again
        //
        // ... and 3 sections in y:
        //    below the box where val = 0
        //    in the box where val = _norm
        //    above the box where val = 0 again
        //
        // Furthermore, we have to calculate the correct values for the pixels on the border.

        int ix_left, ix_right, iy_bottom, iy_top;
        double x_left, x_right, y_bottom, y_top;

        // Find the x edges:
        double tmp = 0.5*width + 0.5;
        ix_left = int(-tmp-x0+1);
        ix_right = int(tmp-x0);

        // If the box goes off the image, it's ok, but it will cause us problems
        // later on if we don't change it.  Just use ix_left = 0.
        if (ix_left < 0) { ix_left = 0; x_left = 1.; }

        // If the whole box is off the image, just zero and return.
        else if (ix_left >= m) { val.setZero(); return; }

        // Normal case: calculate the fractional flux in the edge
        else x_left = tmp+x0+ix_left;

        // Now the right side.
        if (ix_right >= m) { ix_right = m-1; x_right = 1.; }
        else if (ix_right < 0) { val.setZero(); return; }
        else x_right = tmp-x0-ix_right;
        xdbg<<"ix_left = "<<ix_left<<" with partial flux "<<x_left<<std::endl;
        xdbg<<"ix_right = "<<ix_right<<" with partial flux "<<x_right<<std::endl;

        // Repeat for y values
        tmp = 0.5*height + 0.5;
        iy_bottom = int(-tmp-y0+1);
        iy_top = int(tmp-y0);

        if (iy_bottom < 0) { iy_bottom = 0; y_bottom = 1.; }
        else if (iy_bottom >= n) { val.setZero(); return; }
        else y_bottom = tmp+y0+iy_bottom;

        if (iy_top >= n) { iy_top = n-1; y_top = 1.; }
//.........这里部分代码省略.........

void LVector::mBasis(
    const tmv::ConstVectorView<double>& x, const tmv::ConstVectorView<double>& y,
    const tmv::ConstVectorView<double>* invsig,
    tmv::MatrixView<T> psi, int order, double sigma)
{
    assert (y.size()==x.size());
    assert (psi.nrows()==x.size() && psi.ncols()==PQIndex::size(order));

    const int N=order;
    const int npts_full = x.size();

    // It's faster to build the psi matrix in blocks so that more of the matrix stays in
    // L1 cache.  For a (typical) 256 KB L2 cache size, this corresponds to 8 columns in the
    // cache, which is pretty good, since we are usually working on 4 columns at a time,
    // plus either X and Y or 3 Lq vectors.
    const int BLOCKING_FACTOR=4096;

    const int max_npts = std::max(BLOCKING_FACTOR,npts_full);
    tmv::DiagMatrix<double> Rsq_full(max_npts);
    tmv::Matrix<double> A_full(max_npts,2);
    tmv::Matrix<double> tmp_full(max_npts,2);
    tmv::DiagMatrix<double> Lmq_full(max_npts);
    tmv::DiagMatrix<double> Lmqm1_full(max_npts);
    tmv::DiagMatrix<double> Lmqm2_full(max_npts);

    for (int ilo=0; ilo<npts_full; ilo+=BLOCKING_FACTOR) {
        const int ihi = std::min(npts_full, ilo + BLOCKING_FACTOR);
        const int npts = ihi-ilo;

        // Cast arguments as diagonal matrices so we can access
        // vectorized element-by-element multiplication
        tmv::ConstDiagMatrixView<double> X = DiagMatrixViewOf(x.subVector(ilo,ihi));
        tmv::ConstDiagMatrixView<double> Y = DiagMatrixViewOf(y.subVector(ilo,ihi));

        // Get the appropriate portion of our temporary matrices.
        tmv::DiagMatrixView<double> Rsq = Rsq_full.subDiagMatrix(0,npts);
        tmv::MatrixView<double> A = A_full.rowRange(0,npts);
        tmv::MatrixView<double> tmp = tmp_full.rowRange(0,npts);

        // We need rsq values twice, so store them here.
        Rsq = X*X;
        Rsq += Y*Y;

        // This matrix will keep track of real & imag parts
        // of prefactor * exp(-r^2/2) (x+iy)^m / sqrt(m!)

        // Build the Gaussian factor
        for (int i=0; i<npts; i++) A.ref(i,0) = std::exp(-0.5*Rsq(i));
        mBasisHelper<T>::applyPrefactor(A.col(0),sigma);
        A.col(1).setZero();

        // Put 1/sigma factor into every point if doing a design matrix:
        if (invsig) A.col(0) *= tmv::DiagMatrixViewOf(invsig->subVector(ilo,ihi));

        // Assign the m=0 column first:
        psi.col( PQIndex(0,0).rIndex(), ilo,ihi ) = A.col(0);

        // Then ascend m's at q=0:
        for (int m=1; m<=N; m++) {
            int rIndex = PQIndex(m,0).rIndex();
            // Multiply by (X+iY)/sqrt(m), including a factor 2 first time through
            tmp = Y * A;
            A = X * A;
            A.col(0) += tmp.col(1);
            A.col(1) -= tmp.col(0);
            A *= m==1 ? 2. : 1./sqrtn(m);

            psi.subMatrix(ilo,ihi,rIndex,rIndex+2) = mBasisHelper<T>::Asign(m%4) * A;
        }

        // Make three DiagMatrix to hold Lmq's during recurrence calculations
        boost::shared_ptr<tmv::DiagMatrixView<double> > Lmq(
            new tmv::DiagMatrixView<double>(Lmq_full.subDiagMatrix(0,npts)));
        boost::shared_ptr<tmv::DiagMatrixView<double> > Lmqm1(
            new tmv::DiagMatrixView<double>(Lmqm1_full.subDiagMatrix(0,npts)));
        boost::shared_ptr<tmv::DiagMatrixView<double> > Lmqm2(
            new tmv::DiagMatrixView<double>(Lmqm2_full.subDiagMatrix(0,npts)));

        for (int m=0; m<=N; m++) {
            PQIndex pq(m,0);
            int iQ0 = pq.rIndex();
            // Go to q=1:
            pq.incN();
            if (pq.pastOrder(N)) continue;

            {   // q == 1
                const int p = pq.getP();
                const int q = pq.getQ();
                const int iQ = pq.rIndex();

                Lmqm1->setAllTo(1.); // This is Lm0.
                *Lmq = Rsq - (p+q-1.);
                *Lmq *= mBasisHelper<T>::Lsign(1.) / (sqrtn(p)*sqrtn(q));

                if (m==0) {
                    psi.col(iQ,ilo,ihi) = (*Lmq) * psi.col(iQ0,ilo,ihi);
                } else {
                    psi.subMatrix(ilo,ihi,iQ,iQ+2) = (*Lmq) * psi.subMatrix(ilo,ihi,iQ0,iQ0+2);
                }
            }
//.........这里部分代码省略.........