Copy Table2D into TableOld

galsim/__init__.py

@@ -99,7 +99,7 @@
 from .angle import Angle, AngleUnit, _Angle, radians, hours, degrees, arcmin, arcsec
 from .catalog import Catalog, Dict, OutputCatalog
 from .scene import COSMOSCatalog
-from .table import LookupTable, LookupTable2D
+from .table import LookupTable, LookupTable2D, LookupTable2DOld
 
 # Exception and Warning classes
 from .errors import GalSimError, GalSimRangeError, GalSimValueError

galsim/table.py

@@ -643,3 +643,222 @@ def __getstate__(self):
 
     def __setstate__(self, d):
         self.__dict__ = d
+
+
+class LookupTable2DOld(object):
+    def __init__(self, x, y, f, interpolant='linear', edge_mode='raise', constant=0):
+        if edge_mode not in ('raise', 'wrap', 'constant'):
+            raise GalSimValueError("Unknown edge_mode.", edge_mode, ('raise', 'wrap', 'constant'))
+
+        self.x = np.ascontiguousarray(x, dtype=float)
+        self.y = np.ascontiguousarray(y, dtype=float)
+        self.f = np.ascontiguousarray(f, dtype=float)
+
+        dx = np.diff(self.x)
+        dy = np.diff(self.y)
+
+        if not all(dx > 0):
+            raise GalSimValueError("x input grids is not strictly increasing.", x)
+        if not all(dy > 0):
+            raise GalSimValueError("y input grids is not strictly increasing.", y)
+
+        fshape = self.f.shape
+        if fshape != (len(x), len(y)):
+            raise GalSimIncompatibleValuesError(
+                "Shape of f incompatible with lengths of x,y", f=f, x=x, y=y)
+
+        if interpolant not in ('linear', 'ceil', 'floor', 'nearest'):
+            raise GalSimValueError("Unknown interpolant.", interpolant,
+                                   ('linear', 'ceil', 'floor', 'nearest'))
+        self.interpolant = interpolant
+
+
+        self.interpolant = interpolant
+        self.edge_mode = edge_mode
+        self.constant = float(constant)
+
+        if self.edge_mode == 'wrap':
+            # Can wrap if x and y arrays are equally spaced ...
+            if np.allclose(dx, dx[0]) and np.allclose(dy, dy[0]):
+                # Underlying Table2D requires us to extend x, y, and f.
+                self.x = np.append(self.x, self.x[-1]+dx[0])
+                self.y = np.append(self.y, self.y[-1]+dy[0])
+                self.f = np.pad(self.f, [(0,1), (0,1)], mode='wrap')
+            if (all(self.f[0] == self.f[-1]) and all(self.f[:,0] == self.f[:,-1])):
+                self.xperiod = self.x[-1] - self.x[0]
+                self.yperiod = self.y[-1] - self.y[0]
+            else:
+                raise GalSimIncompatibleValuesError(
+                    "Cannot use edge_mode='wrap' unless either x and y are equally "
+                    "spaced or first/last row/column of f are identical.",
+                    edge_mode=edge_mode, x=x, y=y, f=f)
+
+    @lazy_property
+    def _tab(self):
+        with convert_cpp_errors():
+            return _galsim._LookupTable2DOld(self.x.ctypes.data, self.y.ctypes.data,
+                                             self.f.ctypes.data, len(self.x), len(self.y),
+                                             self.interpolant)
+    def getXArgs(self):
+        return self.x
+
+    def getYArgs(self):
+        return self.y
+
+    def getVals(self):
+        return self.f
+
+    def _inbounds(self, x, y):
+        return (np.min(x) >= self.x[0] and np.max(x) <= self.x[-1] and
+                np.min(y) >= self.y[0] and np.max(y) <= self.y[-1])
+
+    def _wrap_args(self, x, y):
+        return ((x-self.x[0]) % self.xperiod + self.x[0],
+                (y-self.y[0]) % self.yperiod + self.y[0])
+
+    @property
+    def _bounds(self):
+        return BoundsD(self.x[0], self.x[-1], self.y[0], self.y[-1])
+
+    def _call_raise(self, x, y):
+        if not self._inbounds(x, y):
+            raise GalSimBoundsError("Extrapolating beyond input range.",
+                                    PositionD(x,y), self._bounds)
+
+        if isinstance(x, numbers.Real):
+            return self._tab.interp(x, y)
+        else:
+            xx = np.ascontiguousarray(x.ravel(), dtype=float)
+            yy = np.ascontiguousarray(y.ravel(), dtype=float)
+            f = np.empty_like(xx, dtype=float)
+            self._tab.interpMany(xx.ctypes.data, yy.ctypes.data, f.ctypes.data, len(xx))
+            f = f.reshape(x.shape)
+            return f
+
+    def _call_wrap(self, x, y):
+        x, y = self._wrap_args(x, y)
+        return self._call_raise(x, y)
+
+    def _call_constant(self, x, y):
+        if isinstance(x, numbers.Real):
+            if self._inbounds(x, y):
+                return self._tab.interp(x, y)
+            else:
+                return self.constant
+        else:
+            x = np.array(x, dtype=float, copy=False)
+            y = np.array(y, dtype=float, copy=False)
+            f = np.empty_like(x, dtype=float)
+            f.fill(self.constant)
+            good = ((x >= self.x[0]) & (x <= self.x[-1]) &
+                    (y >= self.y[0]) & (y <= self.y[-1]))
+            xx = np.ascontiguousarray(x[good].ravel(), dtype=float)
+            yy = np.ascontiguousarray(y[good].ravel(), dtype=float)
+            tmp = np.empty_like(xx, dtype=float)
+            self._tab.interpMany(xx.ctypes.data, yy.ctypes.data, tmp.ctypes.data, len(xx))
+            f[good] = tmp
+            return f
+
+    def __call__(self, x, y):
+        if self.edge_mode == 'raise':
+            return self._call_raise(x, y)
+        elif self.edge_mode == 'wrap':
+            return self._call_wrap(x, y)
+        else: # constant
+            return self._call_constant(x, y)
+
+    def _gradient_raise(self, x, y):
+        if not self._inbounds(x, y):
+            raise GalSimBoundsError("Extrapolating beyond input range.",
+                                    PositionD(x,y), self._bounds)
+
+        if isinstance(x, numbers.Real):
+            grad = np.empty(2, dtype=float)
+            self._tab.gradient(x, y, grad.ctypes.data)
+            return grad[0], grad[1]
+        else:
+            xx = np.ascontiguousarray(x.ravel(), dtype=float)
+            yy = np.ascontiguousarray(y.ravel(), dtype=float)
+            dfdx = np.empty_like(xx)
+            dfdy = np.empty_like(xx)
+            self._tab.gradientMany(xx.ctypes.data, yy.ctypes.data,
+                                   dfdx.ctypes.data, dfdy.ctypes.data, len(xx))
+            dfdx = dfdx.reshape(x.shape)
+            dfdy = dfdy.reshape(x.shape)
+            return dfdx, dfdy
+
+    def _gradient_wrap(self, x, y):
+        x, y = self._wrap_args(x, y)
+        return self._gradient_raise(x, y)
+
+    def _gradient_constant(self, x, y):
+        if isinstance(x, numbers.Real):
+            if self._inbounds(x, y):
+                grad = np.empty(2, dtype=float)
+                self._tab.gradient(float(x), float(y), grad.ctypes.data)
+                return tuple(grad)
+            else:
+                return 0.0, 0.0
+        else:
+            x = np.array(x, dtype=float, copy=False)
+            y = np.array(y, dtype=float, copy=False)
+            dfdx = np.empty_like(x, dtype=float)
+            dfdy = np.empty_like(x, dtype=float)
+            dfdx.fill(0.0)
+            dfdy.fill(0.0)
+            good = ((x >= self.x[0]) & (x <= self.x[-1]) &
+                    (y >= self.y[0]) & (y <= self.y[-1]))
+            xx = np.ascontiguousarray(x[good].ravel(), dtype=float)
+            yy = np.ascontiguousarray(y[good].ravel(), dtype=float)
+            tmp1 = np.empty_like(xx, dtype=float)
+            tmp2 = np.empty_like(xx, dtype=float)
+            self._tab.gradientMany(xx.ctypes.data, yy.ctypes.data,
+                                   tmp1.ctypes.data, tmp2.ctypes.data, len(xx))
+            dfdx[good] = tmp1
+            dfdy[good] = tmp2
+            return dfdx, dfdy
+
+    def gradient(self, x, y):
+        if self.edge_mode == 'raise':
+            return self._gradient_raise(x, y)
+        elif self.edge_mode == 'wrap':
+            return self._gradient_wrap(x, y)
+        else: # constant
+            return self._gradient_constant(x, y)
+
+    def __str__(self):
+        return ("galsim.LookupTable2DOld(x=[%s,...,%s], y=[%s,...,%s], "
+                "f=[[%s,...,%s],...,[%s,...,%s]], interpolant=%r, edge_mode=%r)"%(
+            self.x[0], self.x[-1], self.y[0], self.y[-1],
+            self.f[0,0], self.f[0,-1], self.f[-1,0], self.f[-1,-1],
+            self.edge_mode))
+
+    def __repr__(self):
+        return ("galsim.LookupTable2DOld(x=array(%r), y=array(%r), "
+                "f=array(%r), interpolant=%r, edge_mode=%r, constant=%r)"%(
+            self.x.tolist(), self.y.tolist(), self.f.tolist(), self.interpolant, self.edge_mode,
+            self.constant))
+
+    def __eq__(self, other):
+        return (isinstance(other, LookupTable2DOld) and
+                np.array_equal(self.x,other.x) and
+                np.array_equal(self.y,other.y) and
+                np.array_equal(self.f,other.f) and
+                self.interpolant == other.interpolant and
+                self.edge_mode == other.edge_mode and
+                self.constant == other.constant)
+
+    def __ne__(self, other):
+        return not self.__eq__(other)
+
+    def __hash__(self):
+        return hash(("galsim.LookupTable2DOld", tuple(self.x.ravel()), tuple(self.y.ravel()),
+                    tuple(self.f.ravel()), self.interpolant, self.edge_mode, self.constant))
+
+    def __getstate__(self):
+        d = self.__dict__.copy()
+        d.pop('_tab',None)
+        return d
+
+    def __setstate__(self, d):
+        self.__dict__ = d

include/galsim/TableOld.h

@@ -0,0 +1,69 @@
+/* -*- c++ -*-
+ * Copyright (c) 2012-2018 by the GalSim developers team on GitHub
+ * https://github.com/GalSim-developers
+ *
+ * This file is part of GalSim: The modular galaxy image simulation toolkit.
+ * https://github.com/GalSim-developers/GalSim
+ *
+ * GalSim is free software: redistribution and use in source and binary forms,
+ * with or without modification, are permitted provided that the following
+ * conditions are met:
+ *
+ * 1. Redistributions of source code must retain the above copyright notice, this
+ *    list of conditions, and the disclaimer given in the accompanying LICENSE
+ *    file.
+ * 2. Redistributions in binary form must reproduce the above copyright notice,
+ *    this list of conditions, and the disclaimer given in the documentation
+ *    and/or other materials provided with the distribution.
+ */
+
+#ifndef GalSim_TableOld_H
+#define GalSim_TableOld_H
+
+#include <vector>
+#include <algorithm>
+#include <string>
+#include <sstream>
+#include <stdexcept>
+#include <iostream>
+#include <functional>
+
+#include "Std.h"
+#include "OneDimensionalDeviate.h"
+
+namespace galsim {
+
+    /**
+     * @brief A class to represent lookup tables for a function z = f(x, y).
+     */
+    class Table2DOld
+    {
+    public:
+        enum interpolant { linear, floor, ceil, nearest };
+
+        /// Table from xargs, yargs, vals
+        Table2DOld(const double* xargs, const double* yargs, const double* vals,
+                   int Nx, int Ny, interpolant in);
+
+        /// interp, but exception if beyond bounds
+        double lookup(double x, double y) const;
+
+        /// interp many values at once
+        void interpMany(const double* xvec, const double* yvec, double* valvec, int N) const;
+        void interpManyMesh(const double* xvec, const double* yvec, double* valvec,
+                            int outNx, int outNy) const;
+
+        /// Estimate df/dx, df/dy at a single location
+        void gradient(double x, double y, double& dfdxvec, double& dfdyvec) const;
+
+        /// Estimate many df/dx and df/dy values
+        void gradientMany(const double* xvec, const double* yvec,
+                          double* dfdxvec, double* dfdyvec, int N) const;
+
+    protected:
+        class Table2DOldImpl;
+        shared_ptr<Table2DOldImpl> _pimpl;
+    };
+}
+
+#endif

pysrc/TableOld.cpp

@@ -0,0 +1,75 @@
+/* -*- c++ -*-
+ * Copyright (c) 2012-2018 by the GalSim developers team on GitHub
+ * https://github.com/GalSim-developers
+ *
+ * This file is part of GalSim: The modular galaxy image simulation toolkit.
+ * https://github.com/GalSim-developers/GalSim
+ *
+ * GalSim is free software: redistribution and use in source and binary forms,
+ * with or without modification, are permitted provided that the following
+ * conditions are met:
+ *
+ * 1. Redistributions of source code must retain the above copyright notice, this
+ *    list of conditions, and the disclaimer given in the accompanying LICENSE
+ *    file.
+ * 2. Redistributions in binary form must reproduce the above copyright notice,
+ *    this list of conditions, and the disclaimer given in the documentation
+ *    and/or other materials provided with the distribution.
+ */
+
+#include "PyBind11Helper.h"
+#include "TableOld.h"
+
+namespace galsim {
+
+    static Table2DOld* MakeTable2DOld(size_t ix, size_t iy, size_t ivals, int Nx, int Ny,
+                                const char* interp_c)
+    {
+        const double* x = reinterpret_cast<const double*>(ix);
+        const double* y = reinterpret_cast<const double*>(iy);
+        const double* vals = reinterpret_cast<const double*>(ivals);
+        std::string interp(interp_c);
+
+        Table2DOld::interpolant i = Table2DOld::linear;
+        if (interp == "floor") i = Table2DOld::floor;
+        else if (interp == "ceil") i = Table2DOld::ceil;
+        else if (interp == "nearest") i = Table2DOld::nearest;
+
+        return new Table2DOld(x, y, vals, Nx, Ny, i);
+    }
+
+    static void InterpMany2D(const Table2DOld& Table2DOld, size_t ix, size_t iy, size_t ivals, int N)
+    {
+        const double* x = reinterpret_cast<const double*>(ix);
+        const double* y = reinterpret_cast<const double*>(iy);
+        double* vals = reinterpret_cast<double*>(ivals);
+        Table2DOld.interpMany(x, y, vals, N);
+    }
+
+    static void Gradient(const Table2DOld& Table2DOld, double x, double y, size_t igrad)
+    {
+        double* grad = reinterpret_cast<double*>(igrad);
+        Table2DOld.gradient(x, y, grad[0], grad[1]);
+    }
+
+    static void GradientMany(const Table2DOld& Table2DOld,
+                             size_t ix, size_t iy, size_t idfdx, size_t idfdy, int N)
+    {
+        const double* x = reinterpret_cast<const double*>(ix);
+        const double* y = reinterpret_cast<const double*>(iy);
+        double* dfdx = reinterpret_cast<double*>(idfdx);
+        double* dfdy = reinterpret_cast<double*>(idfdy);
+        Table2DOld.gradientMany(x, y, dfdx, dfdy, N);
+    }
+
+    void pyExportTableOld(PY_MODULE& _galsim)
+    {
+        py::class_<Table2DOld>(GALSIM_COMMA "_LookupTable2DOld" BP_NOINIT)
+            .def(PY_INIT(&MakeTable2DOld))
+            .def("interp", &Table2DOld::lookup)
+            .def("interpMany", &InterpMany2D)
+            .def("gradient", &Gradient)
+            .def("gradientMany", &GradientMany);
+    }
+
+} // namespace galsim

pysrc/files.txt

@@ -27,6 +27,7 @@ Random.cpp
 HSM.cpp
 Integ.cpp
 Table.cpp
+TableOld.cpp
 Interpolant.cpp
 Bessel.cpp
 CDModel.cpp