C++ tmv::MatrixView类代码示例

    void SBMoffat::SBMoffatImpl::fillKValue(tmv::MatrixView<std::complex<double> > val,
                                            double kx0, double dkx, double dkxy,
                                            double ky0, double dky, double dkyx) const
    {
        dbg<<"SBMoffat fillKValue\n";
        dbg<<"kx = "<<kx0<<" + i * "<<dkx<<" + j * "<<dkxy<<std::endl;
        dbg<<"ky = "<<ky0<<" + i * "<<dkyx<<" + j * "<<dky<<std::endl;
        assert(val.stepi() == 1);
        assert(val.canLinearize());
        const int m = val.colsize();
        const int n = val.rowsize();
        typedef tmv::VIt<std::complex<double>,1,tmv::NonConj> It;

        kx0 *= _rD;
        dkx *= _rD;
        dkxy *= _rD;
        ky0 *= _rD;
        dky *= _rD;
        dkyx *= _rD;

        It valit = val.linearView().begin();
        for (int j=0;j<n;++j,kx0+=dkxy,ky0+=dky) {
            double kx = kx0;
            double ky = ky0;
            for (int i=0;i<m;++i,kx+=dkx,ky+=dkyx) {
                double ksq = kx*kx + ky*ky;
                *valit++ = _knorm * (this->*_kV)(ksq);
            }
        }
    }

    void SBGaussian::SBGaussianImpl::fillKValue(tmv::MatrixView<std::complex<double> > val,
                                                double kx0, double dkx, int izero,
                                                double ky0, double dky, int jzero) const
    {
        dbg<<"SBGaussian fillKValue\n";
        dbg<<"kx = "<<kx0<<" + i * "<<dkx<<", izero = "<<izero<<std::endl;
        dbg<<"ky = "<<ky0<<" + j * "<<dky<<", jzero = "<<jzero<<std::endl;
        if (izero != 0 || jzero != 0) {
            xdbg<<"Use Quadrant\n";
            fillKValueQuadrant(val,kx0,dkx,izero,ky0,dky,jzero);
        } else {
            xdbg<<"Non-Quadrant\n";
            assert(val.stepi() == 1);
            const int m = val.colsize();
            const int n = val.rowsize();
            typedef tmv::VIt<double,1,tmv::NonConj> It;

            kx0 *= _sigma;
            dkx *= _sigma;
            ky0 *= _sigma;
            dky *= _sigma;

            tmv::Vector<double> gauss_kx(m);
            It kxit = gauss_kx.begin();
            for (int i=0;i<m;++i,kx0+=dkx) *kxit++ = exp(-0.5 * kx0*kx0);
            tmv::Vector<double> gauss_ky(n);
            It kyit = gauss_ky.begin();
            for (int j=0;j<n;++j,ky0+=dky) *kyit++ = exp(-0.5 * ky0*ky0);

            val = _flux * gauss_kx ^ gauss_ky;
        }
    }

    void SBKolmogorov::SBKolmogorovImpl::fillKValue(tmv::MatrixView<std::complex<double> > val,
                                                    double kx0, double dkx, double dkxy,
                                                    double ky0, double dky, double dkyx) const
    {
        dbg<<"SBKolmogorov fillKValue\n";
        dbg<<"kx = "<<kx0<<" + i * "<<dkx<<" + j * "<<dkxy<<std::endl;
        dbg<<"ky = "<<ky0<<" + i * "<<dkyx<<" + j * "<<dky<<std::endl;
        assert(val.stepi() == 1);
        assert(val.canLinearize());
        const int m = val.colsize();
        const int n = val.rowsize();
        typedef tmv::VIt<std::complex<double>,1,tmv::NonConj> It;

        kx0 *= _inv_k0;
        dkx *= _inv_k0;
        dkxy *= _inv_k0;
        ky0 *= _inv_k0;
        dky *= _inv_k0;
        dkyx *= _inv_k0;

        It valit = val.linearView().begin();
        for (int j=0;j<n;++j,kx0+=dkxy,ky0+=dky) {
            double kx = kx0;
            double ky = ky0;
            for (int i=0;i<m;++i,kx+=dkx,ky+=dkyx) *valit++ = _flux * _info->kValue(kx*kx+ky*ky);
        }
    }

    void SBGaussian::SBGaussianImpl::fillXValue(tmv::MatrixView<double> val,
                                                double x0, double dx, double dxy,
                                                double y0, double dy, double dyx) const
    {
        dbg<<"SBGaussian fillXValue\n";
        dbg<<"x = "<<x0<<" + i * "<<dx<<" + j * "<<dxy<<std::endl;
        dbg<<"y = "<<y0<<" + i * "<<dyx<<" + j * "<<dy<<std::endl;
        assert(val.stepi() == 1);
        assert(val.canLinearize());
        const int m = val.colsize();
        const int n = val.rowsize();
        typedef tmv::VIt<double,1,tmv::NonConj> It;

        x0 *= _inv_sigma;
        dx *= _inv_sigma;
        dxy *= _inv_sigma;
        y0 *= _inv_sigma;
        dy *= _inv_sigma;
        dyx *= _inv_sigma;

        It valit = val.linearView().begin();
        for (int j=0;j<n;++j,x0+=dxy,y0+=dy) {
            double x = x0;
            double y = y0;
            for (int i=0;i<m;++i,x+=dx,y+=dyx)
                *valit++ = _norm * std::exp( -0.5 * (x*x + y*y) );
        }
    }

    void SBMoffat::SBMoffatImpl::fillXValue(tmv::MatrixView<double> val,
                                            double x0, double dx, int izero,
                                            double y0, double dy, int jzero) const
    {
        dbg<<"SBMoffat fillXValue\n";
        dbg<<"x = "<<x0<<" + i * "<<dx<<", izero = "<<izero<<std::endl;
        dbg<<"y = "<<y0<<" + j * "<<dy<<", jzero = "<<jzero<<std::endl;
        if (izero != 0 || jzero != 0) {
            xdbg<<"Use Quadrant\n";
            fillXValueQuadrant(val,x0,dx,izero,y0,dy,jzero);
        } else {
            xdbg<<"Non-Quadrant\n";
            assert(val.stepi() == 1);
            const int m = val.colsize();
            const int n = val.rowsize();
            typedef tmv::VIt<double,1,tmv::NonConj> It;

            x0 *= _inv_rD;
            dx *= _inv_rD;
            y0 *= _inv_rD;
            dy *= _inv_rD;

            for (int j=0;j<n;++j,y0+=dy) {
                double x = x0;
                double ysq = y0*y0;
                It valit = val.col(j).begin();
                for (int i=0;i<m;++i,x+=dx) {
                    double rsq = x*x + ysq;
                    if (rsq > _maxRrD_sq) *valit++ = 0.;
                    else *valit++ = _norm / _pow_beta(1.+rsq, _beta);
                }
            }
        }
    }

    void SBKolmogorov::SBKolmogorovImpl::fillKValue(tmv::MatrixView<std::complex<double> > val,
                                                    double kx0, double dkx, int izero,
                                                    double ky0, double dky, int jzero) const
    {
        dbg<<"SBKolmogorov fillKValue\n";
        dbg<<"kx = "<<kx0<<" + i * "<<dkx<<", izero = "<<izero<<std::endl;
        dbg<<"ky = "<<ky0<<" + j * "<<dky<<", jzero = "<<jzero<<std::endl;
        if (izero != 0 || jzero != 0) {
            xdbg<<"Use Quadrant\n";
            fillKValueQuadrant(val,kx0,dkx,izero,ky0,dky,jzero);
        } else {
            xdbg<<"Non-Quadrant\n";
            assert(val.stepi() == 1);
            const int m = val.colsize();
            const int n = val.rowsize();
            typedef tmv::VIt<std::complex<double>,1,tmv::NonConj> It;

            kx0 *= _inv_k0;
            dkx *= _inv_k0;
            ky0 *= _inv_k0;
            dky *= _inv_k0;

            for (int j=0;j<n;++j,ky0+=dky) {
                double kx = kx0;
                double kysq = ky0*ky0;
                It valit = val.col(j).begin();
                for (int i=0;i<m;++i,kx+=dkx) *valit++ = _flux * _info->kValue(kx*kx + kysq);
            }
        }
    }

    void SBGaussian::SBGaussianImpl::fillXValue(tmv::MatrixView<double> val,
                                                double x0, double dx, int izero,
                                                double y0, double dy, int jzero) const
    {
        dbg<<"SBGaussian fillXValue\n";
        dbg<<"x = "<<x0<<" + i * "<<dx<<", izero = "<<izero<<std::endl;
        dbg<<"y = "<<y0<<" + j * "<<dy<<", jzero = "<<jzero<<std::endl;
        if (izero != 0 || jzero != 0) {
            xdbg<<"Use Quadrant\n";
            fillXValueQuadrant(val,x0,dx,izero,y0,dy,jzero);
        } else {
            xdbg<<"Non-Quadrant\n";
            assert(val.stepi() == 1);
            const int m = val.colsize();
            const int n = val.rowsize();
            typedef tmv::VIt<double,1,tmv::NonConj> It;

            x0 *= _inv_sigma;
            dx *= _inv_sigma;
            y0 *= _inv_sigma;
            dy *= _inv_sigma;

            // The Gaussian profile is separable:
            //    val = _norm * exp(-0.5 * (x*x + y*y)
            //        = _norm * exp(-0.5 * x*x) * exp(-0.5 * y*y)
            tmv::Vector<double> gauss_x(m);
            It xit = gauss_x.begin();
            for (int i=0;i<m;++i,x0+=dx) *xit++ = exp(-0.5 * x0*x0);
            tmv::Vector<double> gauss_y(n);
            It yit = gauss_y.begin();
            for (int j=0;j<n;++j,y0+=dy) *yit++ = exp(-0.5 * y0*y0);

            val = _norm * gauss_x ^ gauss_y;
        }
    }

    void SBKolmogorov::SBKolmogorovImpl::fillXValue(tmv::MatrixView<double> val,
                                                    double x0, double dx, int izero,
                                                    double y0, double dy, int jzero) const
    {
        dbg<<"SBKolmogorov fillXValue\n";
        dbg<<"x = "<<x0<<" + i * "<<dx<<", izero = "<<izero<<std::endl;
        dbg<<"y = "<<y0<<" + j * "<<dy<<", jzero = "<<jzero<<std::endl;
        if (izero != 0 || jzero != 0) {
            xdbg<<"Use Quadrant\n";
            fillXValueQuadrant(val,x0,dx,izero,y0,dy,jzero);
        } else {
            xdbg<<"Non-Quadrant\n";
            assert(val.stepi() == 1);
            const int m = val.colsize();
            const int n = val.rowsize();
            typedef tmv::VIt<double,1,tmv::NonConj> It;

            x0 *= _k0;
            dx *= _k0;
            y0 *= _k0;
            dy *= _k0;

            for (int j=0;j<n;++j,y0+=dy) {
                double x = x0;
                double ysq = y0*y0;
                It valit = val.col(j).begin();
                for (int i=0;i<m;++i,x+=dx) {
                    double r = sqrt(x*x + ysq);
                    *valit++ = _xnorm * _info->xValue(r);
                }
            }
        }
    }

    void SBTopHat::SBTopHatImpl::fillXValue(tmv::MatrixView<double> val,
                                            double x0, double dx, int izero,
                                            double y0, double dy, int jzero) const
    {
        dbg<<"SBTopHat 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;

        val.setZero();
        // The columns to consider have -r0 <= y < r0
        // given that y = y0 + j dy
        double absdx = std::abs(dx);
        double absdy = std::abs(dy);
        int j1 = std::max(0, int(std::ceil(-_r0/absdy - y0/dy)));
        int j2 = std::min(n, int(std::ceil(_r0/absdy - y0/dy)));
        y0 += j1 * dy;
        for (int j=j1;j<j2;++j,y0+=dy) {
            double ysq = y0*y0;
            double xmax = std::sqrt(_r0sq - ysq);
            // Set to _norm all pixels with -xmax <= x < xmax
            // given that x = x0 + i dx.
            int i1 = std::max(0, int(std::ceil(-xmax/absdx - x0/dx)));
            int i2 = std::min(m, int(std::ceil(xmax/absdx - x0/dx)));
            if (i1 < i2)
                val.col(j,i1,i2).setAllTo(_norm);
        }
    }

    void SBBox::SBBoxImpl::fillKValue(tmv::MatrixView<std::complex<double> > val,
                                      double kx0, double dkx, int izero,
                                      double ky0, double dky, int jzero) const
    {
        dbg<<"SBBox fillKValue\n";
        dbg<<"kx = "<<kx0<<" + i * "<<dkx<<", izero = "<<izero<<std::endl;
        dbg<<"ky = "<<ky0<<" + j * "<<dky<<", jzero = "<<jzero<<std::endl;
        if (izero != 0 || jzero != 0) {
            xdbg<<"Use Quadrant\n";
            fillKValueQuadrant(val,kx0,dkx,izero,ky0,dky,jzero);
        } else {
            xdbg<<"Non-Quadrant\n";
            assert(val.stepi() == 1);
            const int m = val.colsize();
            const int n = val.rowsize();
            typedef tmv::VIt<double,1,tmv::NonConj> It;

            kx0 *= _wo2pi;
            dkx *= _wo2pi;
            ky0 *= _ho2pi;
            dky *= _ho2pi;

            // The Box profile in Fourier space is separable:
            //    val(x,y) = _flux * sinc(x * _width/2pi) * sinc(y * _height/2pi)
            tmv::Vector<double> sinc_kx(m);
            It kxit = sinc_kx.begin();
            for (int i=0;i<m;++i,kx0+=dkx) *kxit++ = sinc(kx0);
            tmv::Vector<double> sinc_ky(n);
            It kyit = sinc_ky.begin();
            for (int j=0;j<n;++j,ky0+=dky) *kyit++ = sinc(ky0);

            val = _flux * sinc_kx ^ sinc_ky;
        }
    }

void LVector::kBasis(
    const tmv::ConstVectorView<double>& kx, const tmv::ConstVectorView<double>& ky,
    tmv::MatrixView<std::complex<double> > psi_k, int order, double sigma)
{
    assert(ky.size() == kx.size() && psi_k.nrows() == kx.size());
    assert(psi_k.ncols()==PQIndex::size(order));
    mBasis(kx, ky, 0, psi_k, order, sigma);
}

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

 void SBShapelet::SBShapeletImpl::fillXValue(
     tmv::MatrixView<double> val,
     const tmv::Matrix<double>& x, const tmv::Matrix<double>& y) const
 {
     dbg<<"order = "<<_bvec.getOrder()<<", sigma = "<<_sigma<<std::endl;
     xdbg<<"fillXValue with bvec = "<<_bvec<<std::endl;
     assert(val.stepi() == 1);
     assert(val.canLinearize());
     const int m = val.colsize();
     const int n = val.rowsize();
     tmv::Matrix<double> psi(m*n,_bvec.size());
     LVector::basis(x.constLinearView(),y.constLinearView(),psi.view(),
                    _bvec.getOrder(),_sigma);
     val.linearView() = psi * _bvec.rVector();
 }

 void SBConvolve::SBConvolveImpl::fillKValue(tmv::MatrixView<std::complex<double> > val,
                                             double x0, double dx, double dxy,
                                             double y0, double dy, double dyx) const
 {
     dbg<<"SBConvolve fillKValue\n";
     dbg<<"x = "<<x0<<" + ix * "<<dx<<" + iy * "<<dxy<<std::endl;
     dbg<<"y = "<<y0<<" + ix * "<<dyx<<" + iy * "<<dy<<std::endl;
     ConstIter pptr = _plist.begin();
     assert(pptr != _plist.end());
     GetImpl(*pptr)->fillKValue(val,x0,dx,dxy,y0,dy,dyx);
     if (++pptr != _plist.end()) {
         tmv::Matrix<std::complex<double> > val2(val.colsize(),val.rowsize());
         for (; pptr != _plist.end(); ++pptr) {
             GetImpl(*pptr)->fillKValue(val2.view(),x0,dx,dxy,y0,dy,dyx);
             val = ElemProd(val,val2);
         }
     }
 }

 void SBAdd::SBAddImpl::fillXValue(tmv::MatrixView<double> val,
                                   double x0, double dx, double dxy,
                                   double y0, double dy, double dyx) const
 {
     dbg<<"SBAdd fillXValue\n";
     dbg<<"x = "<<x0<<" + i * "<<dx<<" + j * "<<dxy<<std::endl;
     dbg<<"y = "<<y0<<" + i * "<<dyx<<" + j * "<<dy<<std::endl;
     ConstIter pptr = _plist.begin();
     assert(pptr != _plist.end());
     GetImpl(*pptr)->fillXValue(val,x0,dx,dxy,y0,dy,dyx);
     if (++pptr != _plist.end()) {
         tmv::Matrix<double> val2(val.colsize(),val.rowsize());
         for (; pptr != _plist.end(); ++pptr) {
             GetImpl(*pptr)->fillXValue(val2.view(),x0,dx,dxy,y0,dy,dyx);
             val += val2;
         }
     }
 }