#1265 quinticje

galsim/interpolant.py

@@ -16,7 +16,7 @@
 #    and/or other materials provided with the distribution.
 #
 
-__all__ = [ 'Interpolant', 'Nearest', 'Linear', 'Cubic', 'Quintic',
+__all__ = [ 'Interpolant', 'Nearest', 'Linear', 'Cubic', 'Quintic', 'QuinticBis'
             'Lanczos', 'SincInterpolant', 'Delta', ]
 
 import math
@@ -56,6 +56,7 @@ def from_name(name, tol=None, gsparams=None):
             - 'linear' = `Linear`
             - 'cubic' = `Cubic`
             - 'quintic' = `Quintic`
+            - 'quinticbis' = 'QuinticBis'
             - 'lanczosN' = `Lanczos`  (where N is an integer, given the ``n`` parameter)
 
         In addition, if you want to specify the ``conserve_dc`` option for `Lanczos`, you can
@@ -91,6 +92,8 @@ def from_name(name, tol=None, gsparams=None):
             return Linear(gsparams=gsparams)
         elif name.lower() == 'cubic' :
             return Cubic(gsparams=gsparams)
+        elif name.lower() == 'quinticbis':
+            return QuinticBis(gsparams=gsparams)
         elif name.lower() == 'nearest':
             return Nearest(gsparams=gsparams)
         elif name.lower() == 'delta':
@@ -99,7 +102,7 @@ def from_name(name, tol=None, gsparams=None):
             return SincInterpolant(gsparams=gsparams)
         else:
             raise GalSimValueError("Invalid Interpolant name %s.",name,
-                                   ('linear', 'cubic', 'quintic', 'lanczosN', 'nearest', 'delta',
+                                   ('linear', 'cubic', 'quintic', 'quinticbis', 'lanczosN', 'nearest', 'delta',
                                     'sinc'))
 
     @property
@@ -111,12 +114,12 @@ def gsparams(self):
     @property
     def positive_flux(self):
         """The positive-flux fraction of the interpolation kernel."""
-        return self._i.getPositiveFlux();
+        return self._i.getPositiveFlux()
 
     @property
     def negative_flux(self):
         """The negative-flux fraction of the interpolation kernel."""
-        return self._i.getNegativeFlux();
+        return self._i.getNegativeFlux()
 
     @property
     def tol(self):
@@ -601,6 +604,100 @@ def krange(self):
         return 1.6058208066649935 / self._gsparams.kvalue_accuracy**(1./3.)
 
 
+class QuinticBis(Interpolant):
+    """Fifth order interpolation
+
+    This quintic interpolant is a variation over the 'Quintic' interpolant by Bernstein & Gruen, 
+    found by Campagne Jean-Eric. 
+
+    To summarise the differences between the two quintic kernels, one can say that
+    1. Both kernels share the "quintic" property of iexact interpolation of 
+         fourth-order polynomial functions, and satisfies the partition of unity 
+         both in real ans Fourier spaces.
+    2. :math:`K^{BG}(x)` is more flat than :math:`K^{JE}(x)` around :math `x=0`:
+        .. math::
+            K^{BG}(x) \approx 1. - 7.91667 |x|^3
+        while
+        .. math::
+            K^{JE}(x) \approx 1. - 4.17375 x^2
+    3.  :math:`\widehat{K}^G(u)` is a little more flat than :math:`\widehat{K}^{JE}(u)` around :math:`u=0`:
+        .. math::
+            \widehat{K}^{BG}(u) \approx 1. - 127.168 u^6 
+        while
+        .. math::
+            \widehat{K}^{JE}(u) \approx 1. - 188.312 u^6
+    4. :math:`\widehat{K}^{JE}(u)` is more flat than :math:`\widehat{K}^G(u)` around :math:`$u=1``
+         (ie. :math:`j=1` 1st ghost)
+        .. math::
+            \widehat{K}^{BG}(u) \approx -34.2608 (u-1.)^5
+        while
+        .. math::
+            \widehat{K}^{JE}(u) \approx 70.4229 (u-1.)^6
+    5. :math:`\widehat{K}^{JE}(u)` has higher :math:`j>1` ghosts (ie. coefficient of :math`(u-j)^5`) than 
+    :math:`\widehat{K}^{G}(u)` has a matter of compensation of the 1st ghost intensity. 
+    The ghost intensity of both kernels decrease well as $j$ increase.  
+
+    Parameters:
+        gsparams:   An optional `GSParams` argument. [default: None]
+    """
+    def __init__(self, gsparams=None):
+        self._gsparams = GSParams.check(gsparams)
+
+    @lazy_property
+    def _i(self):
+        return _galsim.QuinticBis(self._gsparams._gsp)
+
+    def __repr__(self):
+        return "galsim.QuinticBis(gsparams=%r)"%(self._gsparams)
+
+    def __str__(self):
+        return "galsim.QuinticBis()"
+
+    @property
+    def xrange(self):
+        """The maximum extent of the interpolant from the origin (in pixels).
+        """
+        return 3.
+
+    @property
+    def ixrange(self):
+        """The total integral range of the interpolant.  Typically 2 * xrange.
+        """
+        return 6
+
+    @property
+    def krange(self):
+        """The maximum extent of the interpolant in Fourier space (in 1/pixels).
+        """
+        # from C++ code:
+        # umax = 1. + std::pow((25.*sqrt(5.)/(135.*pi3-9.*pi5))/tol, 1./3.);
+        # kmax = 2 Pi umax
+        return 2.0*np.pi*(1.0 + 0.33925584826755739773 / self._gsparams.kvalue_accuracy ** (1./3.))
+
+    # numerical values for unit_integrals
+    _unit_integrals = np.array([0.82900606447524252394, 0.10593642386579032927,
+                                -0.023937069546316131708,0.0034976134429045404698,],dtype=float)
+
+    def unit_integrals(self, max_len=None):
+        """Compute the unit integrals of the real-space kernel.
+
+        integrals[i] = int(xval(x), i-0.5, i+0.5)
+
+        Parameters:
+            max_len:    The maximum length of the returned array.  This is usually only relevant
+                        for SincInterpolant, where xrange = inf.
+        Returns:
+            integrals:  An array of unit integrals of length max_len or smaller.
+        """
+        # Using Mathematica:
+        # i0 = int(h(x), x=-0.5..0.5) = 2*int(h(x), x=0..0.5) 
+        #    = (-146715 + 9901 \[Pi]^2)/(11520 (-15 + \[Pi]^2))
+        # i1 = int(h(x), x=-0.5..1.5) = (-62835 + 3829 \[Pi]^2)/(46080 (-15 + \[Pi]^2))
+        # i2 = int(h(x), x=1.5..2.5) = (6195 - 341 \[Pi]^2)/(23040 (-15 + \[Pi]^2))
+        # i3 = int(h(x), x=2.5..3.0) = (-1725 + 91 \[Pi]^2)/(46080 (-15 + \[Pi]^2))
+        return self._unit_integrals
+
+
 class Lanczos(Interpolant):
     """The Lanczos interpolation filter, nominally sinc(x)*sinc(x/n)
 

include/galsim/Interpolant.h

@@ -456,6 +456,53 @@ namespace galsim {
         static std::map<double,double> _cache_umax;
     };
 
+
+
+    /**
+     * @brief Piecewise-quintic polynomial interpolant, ideal for Fourier-space interpolation
+     *
+     * It is a variation on the  Bernstein & Gruen quintic kernel with attetion on the
+     * Taylor expension near u=1
+     */
+
+    class PUBLIC_API QuinticBis : public Interpolant
+    {
+    public:
+        /**
+         * @brief Constructor
+         * @param[in] gsparams GSParams object storing constants that control the accuracy of
+         *                     operations, if different from the default.
+         */
+        QuinticBis(const GSParams& gsparams);
+        ~QuinticBis() {}
+
+        double xrange() const { return _range; }
+        int ixrange() const { return 6; }
+        double urange() const { return _uMax; }
+
+        double xval(double x) const;
+        double uval(double u) const;
+
+        // Override numerical calculation with known analytic integral
+        // Contrary to Bernstein & Gruen kernel, here the roots are located at integer values on real axis
+        double getPositiveFlux() const { return 271./240.; }
+        double getNegativeFlux() const { return 31./240.; }
+
+        std::string makeStr() const;
+
+    private:
+        double _range; // Reduce range slightly from n so we're not using zero-valued endpoints.
+        shared_ptr<TableBuilder> _tab; // Tabulated Fourier transform
+        double _uMax;  // Truncation point for Fourier transform
+
+        // Calculate the FT from a direct integration.
+        double uCalc(double u, double tolerance) const;
+
+        // Store the tables in a map, so repeat constructions are quick.
+        static std::map<double,shared_ptr<TableBuilder> > _cache_tab;
+        static std::map<double,double> _cache_umax;
+    };
+
     /**
      * @brief The Lanczos interpolation filter, nominally sinc(x)*sinc(x/n), truncated at +/-n.
      *

pysrc/Interpolant.cpp

@@ -66,6 +66,9 @@ namespace galsim {
 
         py::class_<Quintic, Interpolant>(_galsim, "Quintic")
             .def(py::init<GSParams>());
+
+        py::class_<QuinticBis, Interpolant>(_galsim, "QuinticBis")
+            .def(py::init<GSParams>());
     }
 
 } // namespace galsim

src/Interpolant.cpp

@@ -556,6 +556,99 @@ namespace galsim {
     }
 
 
+
+    //
+    // QuinticBis
+    //
+
+    double QuinticBis::xval(double x) const
+    {
+        x = std::abs(x);
+        double x2 = x * x;
+        double pi2 = M_PI * M_PI;
+
+        if (x <= 1.)
+            return (15. *(-12. + x2*(27. + x*(-13. + (-3. + x)*x)))
+                + pi2 * (12. - x2*(15. + x*(35. + x*(-63. + 25.*x)))))/(12.*(-15. +pi2));
+        else if (x <= 2.)
+            return ((-2. + x) * (-1 + x)*(-15.*(24. + x*(-3. + (-6. + x)*x))
+                + pi2 * (-48. + x * (153. + x*(-114. + 25.*x)))))/(24.*(-15. + pi2));
+        else if (x <= 3.)
+            return -(((-3. + x)*(-3.+x)*(-2. + x)*(-3.* (-7. + x)* x
+                + pi2*(-3. + x)*(-8. + 5.*x)))/(24.*(-15. + pi2)));
+        else
+            return 0.;
+    }
+
+    double QuinticBis::uval(double u) const
+    {
+        u = std::abs(u);
+#ifdef USE_TABLES
+        return u>_uMax ? 0. : (*_tab)(u);
+#else
+        double pi2 = M_PI * M_PI;
+        double piu = M_PI * u;
+        double scu = math::sinc(u);
+        double cu = cos(piu);
+        double su = sin(piu);
+        double ssq = scu * scu;
+        return (scu*ssq*ssq * (M_PI*(24.*M_PI* (-1. + u*u)*cu - (39. + 7.*pi2)* u*su)
+                + 5.*(-3. + 5.*pi2)*scu))/(-15. + pi2);
+#endif
+    }
+
+    class QuinticBisIntegrand
+    {
+    public:
+        QuinticBisIntegrand(double u, const QuinticBis& q): _u(u), _q(q) {}
+        double operator()(double x) const { return _q.xval(x)*std::cos(2*M_PI*_u*x); }
+    private:
+        double _u;
+        const QuinticBis& _q;
+    };
+
+    double QuinticBis::uCalc(double u, double tolerance) const
+    {
+        QuinticBisIntegrand qi(u, *this);
+        return 2.*( integ::int1d(qi, 0., 1., 0.1*tolerance, 0.1*tolerance)
+                    + integ::int1d(qi, 1., 2., 0.1*tolerance, 0.1*tolerance)
+                    + integ::int1d(qi, 2., 3., 0.1*tolerance, 0.1*tolerance));
+    }
+
+    QuinticBis::QuinticBis(const GSParams& gsparams) : Interpolant(gsparams)
+    {
+        dbg<<"Start QuinticBis\n";
+        _range = 3.;
+        double pi2 = M_PI * M_PI;
+        double pi3 = M_PI * pi2;
+        double pi5 = pi2 * pi3;
+        // uMax is the value where |ft| <= tolerance
+        // ft = (Sinc[u]^5 (24 pi^2 (1 - u^2) Cos[pi u] 
+        //               + (15 + pi^2 (-25 + (39 + 7 pi^2) u^2)) Sinc[pi u]))/(15 - pi^2)
+        // Mathematica: It turns out that the function  
+        // g(u) = 24 pi^2 /((15 - pi^2) pi^5) (25 Sqrt[5])/216 Abs[1/ (u - 1)^3]
+        // gives a good approximation of max(|ft|) on each [i,i+1] intervalle, i integer
+        // nb. Max[Sinc[u]^5 Cos[pi u]] = (25 Sqrt[5])/216
+        _uMax = 1. + std::pow((25.*sqrt(5.)/(135.*pi3-9.*pi5))/gsparams.kvalue_accuracy, 1./3.);
+    }
+
+    std::map<double,shared_ptr<TableBuilder> > QuinticBis::_cache_tab;
+    std::map<double,double> QuinticBis::_cache_umax;
+
+    // nb. Contrary to Quintic from Bernstein & Gruen, 
+    // for QuinticBis the change of sign occurs on x=i with i non-zero integer
+    // So the default sampler is ok.
+
+    std::string QuinticBis::makeStr() const
+    {
+        std::ostringstream oss(" ");
+        oss.precision(std::numeric_limits<double>::digits10 + 4);
+        oss << "galsim._galsim.QuinticBis(";
+        oss << "galsim._galsim.GSParams("<<_gsparams<<"))";
+        return oss.str();
+    }
+
+
     //
     // Lanczos
     //

tests/interpolant_comparison_files/absfKQuinticBis_test.txt

@@ -0,0 +1,3 @@
+8.465317820709420e-09
+8.877555895048230e-01
+1.864649926651702e-09