From 3cdd034a6a477196d134812bcc2aff85d1c6b6e0 Mon Sep 17 00:00:00 2001 From: inakleinbottle <41870650+inakleinbottle@users.noreply.github.com> Date: Fri, 24 Nov 2023 13:43:02 +0000 Subject: [PATCH] Fix conversion to numpy (#52) * Added tests for to_array functionality * Stream out operator for TypeInfo * formatting * basis function for algebra types * enhanced framework for converting from algebra to array. * remove some redundant asserts * remove cvref on type deductions * formatting * always import numpy if available * lots of little details * moved some functions around to make them accessible * formatting * formatting * imort numpy function * moved definition of __array__ to numpy files * formatting * Create zeroed array first * Added missing zeros to expected arrays * no ownership test - everything copies for now * Added test for array of rationals * fixed some issues with scalar allocations * Updated changelog * Added test for poly scalars to_array * Fixed missing break in switch statement. --- CHANGELOG | 2 + .../include/roughpy/algebra/algebra_base.h | 3 +- .../roughpy/algebra/algebra_base_impl.h | 9 + device/CMakeLists.txt | 1 + device/include/roughpy/device/core.h | 2 +- device/src/core.cpp | 41 +++ roughpy/src/algebra/setup_algebra_type.h | 16 +- roughpy/src/args/numpy.cpp | 309 ++++++++++++++---- roughpy/src/args/numpy.h | 57 +++- roughpy/src/roughpy_module.cpp | 5 + roughpy/src/scalars/r_py_polynomial.cpp | 3 - roughpy/src/scalars/r_py_polynomial.h | 4 + roughpy/src/scalars/scalars.cpp | 2 +- roughpy/src/scalars/scalars.h | 2 + scalars/include/roughpy/scalars/scalar.h | 8 +- scalars/src/scalar.cpp | 9 +- .../src/types/aprational/ap_rational_type.cpp | 11 +- .../src/types/aprational/ap_rational_type.h | 4 + tests/algebra/test_free_tensor.py | 90 ++++- 19 files changed, 470 insertions(+), 108 deletions(-) create mode 100644 device/src/core.cpp diff --git a/CHANGELOG b/CHANGELOG index 2feeb85c..e7c4e5ac 100644 --- a/CHANGELOG +++ b/CHANGELOG @@ -1,6 +1,8 @@ Version 0.1.0: - Added framework for integrating device support and redesigned scalars module to accommodate the changes. - Made changes to type deduction in constructors to avoid exceptions when providing lists of python ints/floats. + - Changed the implementation of __array__ for algebra types. A copy is always made, and the size of the array is always + equal to the dimension of the chosen width/depth composition. Version 0.0.8: diff --git a/algebra/include/roughpy/algebra/algebra_base.h b/algebra/include/roughpy/algebra/algebra_base.h index b7aee4e2..cbb3f00a 100644 --- a/algebra/include/roughpy/algebra/algebra_base.h +++ b/algebra/include/roughpy/algebra/algebra_base.h @@ -191,7 +191,8 @@ class AlgebraBase return p_impl.get(); } - + RPY_NO_DISCARD + basis_type basis() const; RPY_NO_DISCARD dimn_t dimension() const; RPY_NO_DISCARD diff --git a/algebra/include/roughpy/algebra/algebra_base_impl.h b/algebra/include/roughpy/algebra/algebra_base_impl.h index aea35466..0c385af0 100644 --- a/algebra/include/roughpy/algebra/algebra_base_impl.h +++ b/algebra/include/roughpy/algebra/algebra_base_impl.h @@ -254,6 +254,15 @@ context_pointer AlgebraBase::context() const noexcept return (p_impl) ? p_impl->context() : nullptr; } + +template < + typename Interface, + template class> class DerivedImpl> +typename AlgebraBase::basis_type AlgebraBase::basis() const +{ + return (p_impl) ? p_impl->basis() : basis_type(); +} + template < typename Interface, template class> class DerivedImpl> diff --git a/device/CMakeLists.txt b/device/CMakeLists.txt index ed38b82e..ee5366b0 100644 --- a/device/CMakeLists.txt +++ b/device/CMakeLists.txt @@ -71,6 +71,7 @@ add_roughpy_component(Device src/opencl/ocl_version.h src/buffer.cpp src/buffer_interface.cpp + src/core.cpp src/device_handle.cpp src/device_interface_base.cpp src/device_provider.cpp diff --git a/device/include/roughpy/device/core.h b/device/include/roughpy/device/core.h index 16c17dbc..ad034230 100644 --- a/device/include/roughpy/device/core.h +++ b/device/include/roughpy/device/core.h @@ -403,7 +403,7 @@ constexpr TypeInfo type_info() noexcept 1U}; } - +std::ostream& operator<<(std::ostream& os, const TypeInfo& code); RPY_SERIAL_SERIALIZE_FN_EXT(TypeInfo) diff --git a/device/src/core.cpp b/device/src/core.cpp new file mode 100644 index 00000000..665a146d --- /dev/null +++ b/device/src/core.cpp @@ -0,0 +1,41 @@ +// +// Created by sam on 11/23/23. +// +#include "core.h" + + +std::ostream& rpy::devices::operator<<(std::ostream& os, const TypeInfo& info) +{ + switch (info.code) { + case TypeCode::Int: os << "int" << info.bytes * CHAR_BIT; + break; + case TypeCode::UInt: os << "uint" << info.bytes * CHAR_BIT; + break; + case TypeCode::Float: os << "float" << info.bytes * CHAR_BIT; + break; + case TypeCode::OpaqueHandle: os << "opaque"; + break; + case TypeCode::BFloat: os << "bfloat" << info.bytes * CHAR_BIT; + break; + case TypeCode::Complex: os << "complex" << info.bytes * CHAR_BIT; + break; + case TypeCode::Bool: os << "bool"; + break; + case TypeCode::Rational: os << "Rational"; + break; + case TypeCode::ArbitraryPrecisionInt: os << "mp_int"; + break; + case TypeCode::ArbitraryPrecisionUInt: os << "mp_uint"; + break; + case TypeCode::ArbitraryPrecisionFloat: os << "mp_float"; + break; + case TypeCode::ArbitraryPrecisionComplex: os << "mp_complex"; + break; + case TypeCode::ArbitraryPrecisionRational: os << "Rational"; + break; + case TypeCode::APRationalPolynomial: os << "PolyRational"; + break; + } + + return os; +} diff --git a/roughpy/src/algebra/setup_algebra_type.h b/roughpy/src/algebra/setup_algebra_type.h index a580a057..b63135a8 100644 --- a/roughpy/src/algebra/setup_algebra_type.h +++ b/roughpy/src/algebra/setup_algebra_type.h @@ -225,20 +225,8 @@ void setup_algebra_type(py::class_& klass) // setup conversion to numpy array #ifdef ROUGHPY_WITH_NUMPY klass.def("__array__", [](const Alg& self) { - // py::dtype dtype = dtype_from(self.coeff_type()); - py::dtype dtype = ctype_to_npy_dtype(self.coeff_type()); - - auto dense_data = self.dense_data(); - if (dense_data) { - const auto dense_data_inner = *dense_data; - return py::array( - dtype, - {dense_data_inner.size()}, - {}, - dense_data_inner.pointer() - ); - } - return py::array(dtype); + return algebra_to_array(self); + // return py::array(); }); #endif diff --git a/roughpy/src/args/numpy.cpp b/roughpy/src/args/numpy.cpp index 93b63b97..6a42509d 100644 --- a/roughpy/src/args/numpy.cpp +++ b/roughpy/src/args/numpy.cpp @@ -27,62 +27,119 @@ // POSSIBILITY OF SUCH DAMAGE. #include "numpy.h" +#include #include #include +#include + +#define NPY_NO_DEPRECATED_API NPY_1_17_API_VERSION +#include + +#include "scalars/r_py_polynomial.h" +#include "scalars/scalars.h" + +#include using namespace rpy; -namespace { +// namespace { +// +// enum NpyDtypes : int +// { +// NPY_BOOL = 0, +// NPY_BYTE, +// NPY_UBYTE, +// NPY_SHORT, +// NPY_USHORT, +// NPY_INT, +// NPY_UINT, +// NPY_LONG, +// NPY_ULONG, +// NPY_LONGLONG, +// NPY_ULONGLONG, +// NPY_FLOAT, +// NPY_DOUBLE, +// NPY_LONGDOUBLE, +// NPY_CFLOAT, +// NPY_CDOUBLE, +// NPY_CLONGDOUBLE, +// NPY_OBJECT = 17, +// NPY_STRING, +// NPY_UNICODE, +// NPY_VOID, +// /* +// * New 1.6 types appended, may be integrated +// * into the above in 2.0. +// */ +// NPY_DATETIME, +// NPY_TIMEDELTA, +// NPY_HALF, +// +// NPY_NTYPES, +// NPY_NOTYPE, +// NPY_CHAR, +// NPY_USERDEF = 256, /* leave room for characters */ +// +// /* The number of types not including the new 1.6 types */ +// NPY_NTYPES_ABI_COMPATIBLE = 21 +// }; +// +// }// namespace -enum NpyDtypes : int +static inline int info_to_typenum(const devices::TypeInfo& info) { - NPY_BOOL = 0, - NPY_BYTE, - NPY_UBYTE, - NPY_SHORT, - NPY_USHORT, - NPY_INT, - NPY_UINT, - NPY_LONG, - NPY_ULONG, - NPY_LONGLONG, - NPY_ULONGLONG, - NPY_FLOAT, - NPY_DOUBLE, - NPY_LONGDOUBLE, - NPY_CFLOAT, - NPY_CDOUBLE, - NPY_CLONGDOUBLE, - NPY_OBJECT = 17, - NPY_STRING, - NPY_UNICODE, - NPY_VOID, - /* - * New 1.6 types appended, may be integrated - * into the above in 2.0. - */ - NPY_DATETIME, - NPY_TIMEDELTA, - NPY_HALF, - - NPY_NTYPES, - NPY_NOTYPE, - NPY_CHAR, - NPY_USERDEF = 256, /* leave room for characters */ - - /* The number of types not including the new 1.6 types */ - NPY_NTYPES_ABI_COMPATIBLE = 21 -}; - -}// namespace + switch (info.code) { + case devices::TypeCode::Int: switch (info.bytes) { + case 1: return NPY_INT8; + case 2: return NPY_INT16; + case 4: return NPY_INT32; + case 8: return NPY_INT64; + default: break; + } + break; + case devices::TypeCode::UInt: switch (info.bytes) { + case 1: return NPY_UINT8; + case 2: return NPY_UINT16; + case 4: return NPY_UINT32; + case 8: return NPY_UINT64; + default: break; + } + break; + case devices::TypeCode::Float: switch (info.bytes) { + case 2: return NPY_FLOAT16; + case 4: return NPY_FLOAT32; + case 8: return NPY_FLOAT64; + default: break; + } + break; + case devices::TypeCode::OpaqueHandle: break; + case devices::TypeCode::BFloat: break; + case devices::TypeCode::Complex: break; + case devices::TypeCode::Bool: break; + case devices::TypeCode::Rational: + case devices::TypeCode::ArbitraryPrecisionInt: + case devices::TypeCode::ArbitraryPrecisionUInt: + case devices::TypeCode::ArbitraryPrecisionFloat: + case devices::TypeCode::ArbitraryPrecisionComplex: + case devices::TypeCode::ArbitraryPrecisionRational: + case devices::TypeCode::APRationalPolynomial: return NPY_OBJECT; + } + + std::stringstream ss; + ss << "scalar type " << info << " is not supported in conversions"; + + RPY_THROW(std::runtime_error, ss.str()); +} const scalars::ScalarType* python::npy_dtype_to_ctype(pybind11::dtype dtype) { const scalars::ScalarType* type = nullptr; switch (dtype.num()) { - case NPY_FLOAT: type = *scalars::ScalarType::of(); break; - case NPY_DOUBLE: type = *scalars::ScalarType::of(); break; + case NPY_FLOAT: type = *scalars::ScalarType::of(); + break; + case NPY_DOUBLE: type = *scalars::ScalarType::of(); + break; default: // Default behaviour, promote to double type = *scalars::ScalarType::of(); @@ -91,12 +148,10 @@ const scalars::ScalarType* python::npy_dtype_to_ctype(pybind11::dtype dtype) return type; } + pybind11::dtype python::ctype_to_npy_dtype(const scalars::ScalarType* type) { - if (type == *scalars::ScalarType::of()) { return py::dtype("d"); } - if (type == *scalars::ScalarType::of()) { return py::dtype("f"); } - - RPY_THROW(py::type_error, "unsupported data type"); + return py::dtype(info_to_typenum(type->type_info())); } string python::npy_dtype_to_identifier(pybind11::dtype dtype) @@ -104,35 +159,159 @@ string python::npy_dtype_to_identifier(pybind11::dtype dtype) string identifier; switch (dtype.num()) { - case NPY_FLOAT: identifier = "f32"; break; - case NPY_DOUBLE: identifier = "f64"; break; - case NPY_INT: identifier = "i32"; break; - case NPY_UINT: identifier = "u32"; break; + case NPY_FLOAT: identifier = "f32"; + break; + case NPY_DOUBLE: identifier = "f64"; + break; + case NPY_INT: identifier = "i32"; + break; + case NPY_UINT: identifier = "u32"; + break; case NPY_LONG: { - if (sizeof(long) == sizeof(int)) { - identifier = "i32"; - } else { + if (sizeof(long) == sizeof(int)) { identifier = "i32"; } else { identifier = "i64"; } - } break; + } + break; case NPY_ULONG: { - if (sizeof(long) == sizeof(int)) { - identifier = "ui32"; - } else { + if (sizeof(long) == sizeof(int)) { identifier = "ui32"; } else { identifier = "u64"; } - } break; - case NPY_LONGLONG: identifier = "i64"; break; - case NPY_ULONGLONG: identifier = "u64"; break; + } + break; + case NPY_LONGLONG: identifier = "i64"; + break; + case NPY_ULONGLONG: identifier = "u64"; + break; case NPY_BOOL: - case NPY_BYTE: identifier = "i8"; break; - case NPY_UBYTE: identifier = "u8"; break; - case NPY_SHORT: identifier = "i16"; break; - case NPY_USHORT: identifier = "u16"; break; + case NPY_BYTE: identifier = "i8"; + break; + case NPY_UBYTE: identifier = "u8"; + break; + case NPY_SHORT: identifier = "i16"; + break; + case NPY_USHORT: identifier = "u16"; + break; default: RPY_THROW(py::type_error, "unsupported dtype"); } return identifier; } + +namespace { + +template +void write_type_as_py_object(PyObject** dst, Slice data) +{ + const auto type_o = scalars::scalar_type_of(); + RPY_CHECK(type_o); + for (dimn_t i = 0; i < data.size(); ++i) { + Py_XDECREF(dst[i]); + dst[i] = py::cast(scalars::Scalar(*type_o, data[i])).release().ptr(); + } +} + +void write_type_as_py_object(PyObject** dst, + Slice data) +{ + for (dimn_t i = 0; i < data.size(); ++i) { + Py_XDECREF(dst[i]); + dst[i] = PyPolynomial_FromPolynomial( + scalars::rational_poly_scalar(data[i])); + } +} + + +} + + +py::array python::dtl::dense_data_to_array(const scalars::ScalarArray& data, + dimn_t dimension) +{ + RPY_DBG_ASSERT(data.size() <= dimension); + const auto type_info = data.type_info(); + + py::dtype dtype(info_to_typenum(type_info)); + + auto result = python::dtl::new_zero_array_for_stype(*data.type(), dimension); + + if (scalars::traits::is_fundamental(type_info)) { + scalars::dtl::scalar_convert_copy(result.mutable_data(), + type_info, + data.pointer(), + type_info, + data.size()); + } else { + PyArray_FillObjectArray(reinterpret_cast(result.ptr()), + Py_None); + + auto** raw = static_cast(result.mutable_data()); + + switch (type_info.code) { + case devices::TypeCode::ArbitraryPrecisionRational + : write_type_as_py_object( + raw, + data.template as_slice()); + break; + case devices::TypeCode::APRationalPolynomial + : write_type_as_py_object(raw, + data.template as_slice< + devices::rational_poly_scalar>()); + break; + case devices::TypeCode::Int: + case devices::TypeCode::UInt: + case devices::TypeCode::Float: + case devices::TypeCode::OpaqueHandle: + case devices::TypeCode::BFloat: + case devices::TypeCode::Complex: + case devices::TypeCode::Bool: + case devices::TypeCode::Rational: + case devices::TypeCode::ArbitraryPrecisionInt: + case devices::TypeCode::ArbitraryPrecisionUInt: + case devices::TypeCode::ArbitraryPrecisionFloat: + case devices::TypeCode::ArbitraryPrecisionComplex: + RPY_THROW(std::runtime_error, + "this case should have been handled elsewhere"); + } + } + + return result; +} + +py::array python::dtl::new_zero_array_for_stype(const scalars::ScalarType* type, + dimn_t dimension) +{ + auto typenum = info_to_typenum(type->type_info()); + py::dtype dtype(typenum); + py::array result(dtype, {static_cast(dimension)}, {}); + + if (typenum == NPY_OBJECT) { + PyArray_FillObjectArray(reinterpret_cast(result.ptr()), + Py_None); + } else { + PyArray_FILLWBYTE(reinterpret_cast(result.ptr()), 0); + } + + return result; +} + + +void python::dtl::write_entry_to_array(py::array& array, + dimn_t index, + const scalars::Scalar& arg) {} + + +static inline PyTypeObject* import_numpy_impl() noexcept +{ + import_array(); + return &PyArray_Type; +} + +void python::import_numpy() +{ +#ifdef ROUGHPY_WITH_NUMPY + if (import_numpy_impl() == nullptr) { throw py::error_already_set(); } +#endif +} diff --git a/roughpy/src/args/numpy.h b/roughpy/src/args/numpy.h index 0e341edd..0dd28d6f 100644 --- a/roughpy/src/args/numpy.h +++ b/roughpy/src/args/numpy.h @@ -33,7 +33,8 @@ # include "roughpy_module.h" # include -# include +# include +#include # include @@ -46,8 +47,62 @@ py::dtype ctype_to_npy_dtype(const scalars::ScalarType* type); string npy_dtype_to_identifier(py::dtype dtype); +namespace dtl { + +RPY_NO_DISCARD +py::array dense_data_to_array(const scalars::ScalarArray& data, + dimn_t dimension); +RPY_NO_DISCARD +py::array new_zero_array_for_stype(const scalars::ScalarType* type, + dimn_t dimension); +void write_entry_to_array(py::array& array, + dimn_t index, + const scalars::Scalar& arg); + + +} + +template class> + class DerivedImpl> +RPY_NO_DISCARD +inline py::array algebra_to_array( + const algebra::AlgebraBase& alg) +{ + const auto* stype = alg.coeff_type(); + const auto basis = alg.basis(); + const auto dimension = basis.dimension(); + + auto dense_data = alg.dense_data(); + + if (dense_data && dense_data->size() == dimension) { + // Dense and full dimension, borrow + auto dtype = ctype_to_npy_dtype(stype); + return py::array(dtype, {dimension}, {}, dense_data->pointer()); + } + if (dense_data) { return dtl::dense_data_to_array(*dense_data, dimension); } + + auto result = dtl::new_zero_array_for_stype(stype, dimension); + + for (auto&& item : alg) { + dimn_t index = basis.key_to_index(item.key()); + dtl::write_entry_to_array(result, index, item.value()); + } + + return result; +} + }// namespace python }// namespace rpy #endif// ROUGHPY_WITH_NUMPY + +namespace rpy { +namespace python { + +void import_numpy(); + +} +} + + #endif// RPY_PY_ARGS_NUMPY_H_ diff --git a/roughpy/src/roughpy_module.cpp b/roughpy/src/roughpy_module.cpp index 2a399954..2b3d3337 100644 --- a/roughpy/src/roughpy_module.cpp +++ b/roughpy/src/roughpy_module.cpp @@ -43,6 +43,8 @@ #include "scalars/scalars.h" #include "streams/streams.h" +#include "args/numpy.h" + #ifndef ROUGHPY_VERSION_STRING # define ROUGHPY_VERSION_STRING "1.0.0" @@ -63,4 +65,7 @@ PYBIND11_MODULE(_roughpy, m) init_algebra(m); init_streams(m); // init_recombine(m); + + + import_numpy(); } diff --git a/roughpy/src/scalars/r_py_polynomial.cpp b/roughpy/src/scalars/r_py_polynomial.cpp index d6e5b64a..3102c1a6 100644 --- a/roughpy/src/scalars/r_py_polynomial.cpp +++ b/roughpy/src/scalars/r_py_polynomial.cpp @@ -369,9 +369,6 @@ static PyObject* polynomial___setstate__(PyObject* self, PyObject* args) } } -PyObject* PyPolynomial_FromPolynomial(scalars::rational_poly_scalar&& poly -) noexcept; - static PyMethodDef RPyPolynomial_methods[] = { { "degree",(PyCFunction) &polynomial_degree,METH_NOARGS, nullptr }, {"__getstate__", diff --git a/roughpy/src/scalars/r_py_polynomial.h b/roughpy/src/scalars/r_py_polynomial.h index 35fb3398..14a70b7e 100644 --- a/roughpy/src/scalars/r_py_polynomial.h +++ b/roughpy/src/scalars/r_py_polynomial.h @@ -54,6 +54,10 @@ struct RPyPolynomial { extern PyTypeObject RPyPolynomial_Type; +PyObject* PyPolynomial_FromPolynomial(rpy::scalars::rational_poly_scalar&& poly +) noexcept; + + inline bool RPyMonomial_Check(PyObject* obj) { return Py_TYPE(obj) == &RPyMonomial_Type; diff --git a/roughpy/src/scalars/scalars.cpp b/roughpy/src/scalars/scalars.cpp index 06d3f9a6..d34ee3e8 100644 --- a/roughpy/src/scalars/scalars.cpp +++ b/roughpy/src/scalars/scalars.cpp @@ -653,7 +653,7 @@ scalars::KeyScalarArray python::py_to_buffer( } } } - } else { + } else if (object.is_none()) {} else { RPY_THROW( std::invalid_argument, "could not parse argument to a valid scalar array type" diff --git a/roughpy/src/scalars/scalars.h b/roughpy/src/scalars/scalars.h index 569f87c0..633fba52 100644 --- a/roughpy/src/scalars/scalars.h +++ b/roughpy/src/scalars/scalars.h @@ -134,6 +134,8 @@ ArgSizeInfo compute_size_and_type( py::handle arg ); + + void init_scalars(py::module_& m); }// namespace python diff --git a/scalars/include/roughpy/scalars/scalar.h b/scalars/include/roughpy/scalars/scalar.h index 24072c91..0e128b87 100644 --- a/scalars/include/roughpy/scalars/scalar.h +++ b/scalars/include/roughpy/scalars/scalar.h @@ -210,7 +210,7 @@ class RPY_EXPORT Scalar trivial_bytes, type_info(), &value, - devices::type_info() + devices::type_info>() ); } @@ -229,9 +229,11 @@ class RPY_EXPORT Scalar { allocate_data(); auto this_info = type_info(); - auto value_info = devices::type_info(); + auto value_info = devices::type_info>(); if (this_info == value_info) { - construct_inplace(opaque_pointer, std::forward(value)); + construct_inplace( + static_cast*>(opaque_pointer), + std::forward(value)); } else { dtl::scalar_convert_copy( opaque_pointer, diff --git a/scalars/src/scalar.cpp b/scalars/src/scalar.cpp index 404aa3dd..284bbdab 100644 --- a/scalars/src/scalar.cpp +++ b/scalars/src/scalar.cpp @@ -184,8 +184,12 @@ Scalar::Scalar(Scalar&& other) noexcept break; case dtl::ScalarContentType::OpaquePointer: case dtl::ScalarContentType::ConstOpaquePointer: + opaque_pointer = other.opaque_pointer; + other.opaque_pointer = nullptr; + break; case dtl::ScalarContentType::OwnedPointer: opaque_pointer = other.opaque_pointer; + other.opaque_pointer = nullptr; break; case dtl::ScalarContentType::Interface: case dtl::ScalarContentType::OwnedInterface: @@ -198,7 +202,10 @@ Scalar::~Scalar() switch (p_type_and_content_type.get_enumeration()) { case dtl::ScalarContentType::OwnedPointer: RPY_DBG_ASSERT(p_type_and_content_type.is_pointer()); - p_type_and_content_type->free_single(opaque_pointer); + if (opaque_pointer != nullptr) { + p_type_and_content_type->free_single(opaque_pointer); + opaque_pointer = nullptr; + } break; case dtl::ScalarContentType::Interface: case dtl::ScalarContentType::OwnedInterface: interface.~unique_ptr(); diff --git a/scalars/src/types/aprational/ap_rational_type.cpp b/scalars/src/types/aprational/ap_rational_type.cpp index 9db04533..7d2879fb 100644 --- a/scalars/src/types/aprational/ap_rational_type.cpp +++ b/scalars/src/types/aprational/ap_rational_type.cpp @@ -30,11 +30,18 @@ ScalarArray APRationalType::allocate(dimn_t count) const } void* APRationalType::allocate_single() const { - return ScalarType::allocate_single(); + guard_type access(m_lock); + auto [pos, inserted] = m_allocations.insert(new rational_scalar_type()); + RPY_DBG_ASSERT(inserted); + return *pos; } void APRationalType::free_single(void* ptr) const { - ScalarType::free_single(ptr); + guard_type access(m_lock); + auto found = m_allocations.find(ptr); + RPY_CHECK(found != m_allocations.end()); + delete static_cast(ptr); + m_allocations.erase(found); } void APRationalType::convert_copy(ScalarArray& dst, const ScalarArray& src) const diff --git a/scalars/src/types/aprational/ap_rational_type.h b/scalars/src/types/aprational/ap_rational_type.h index 00a3766a..137d8289 100644 --- a/scalars/src/types/aprational/ap_rational_type.h +++ b/scalars/src/types/aprational/ap_rational_type.h @@ -8,11 +8,15 @@ #include "scalar_type.h" #include "scalar_types.h" +#include + namespace rpy { namespace scalars { class APRationalType : public ScalarType { + mutable std::unordered_set m_allocations; + public: APRationalType(); diff --git a/tests/algebra/test_free_tensor.py b/tests/algebra/test_free_tensor.py index 3745cc73..0c892c89 100644 --- a/tests/algebra/test_free_tensor.py +++ b/tests/algebra/test_free_tensor.py @@ -45,8 +45,7 @@ def test_create_single_float_no_ctype(): result = FreeTensor(1.0, width=2, depth=2) np_result = np.array(result) - assert np_result.shape == (1,) - assert_array_equal(np_result, np.array([1.])) + assert_array_equal(np_result, np.array([1., 0, 0, 0, 0, 0, 0])) assert result.width == 2 assert result.max_degree == 2 @@ -55,8 +54,8 @@ def test_create_single_int_no_ctype(): result = FreeTensor(1, width=2, depth=2) np_result = np.array(result) - assert np_result.shape == (1,) - assert_array_equal(np_result, np.array([1.])) + assert np_result.flags.owndata + assert_array_equal(np_result, np.array([1., 0, 0, 0, 0, 0, 0])) assert result.width == 2 assert result.max_degree == 2 @@ -65,9 +64,9 @@ def test_create_list_floats_no_ctype(): result = FreeTensor([0., 1., 2., 3.], width=3, depth=2) np_result = np.array(result) - assert np_result.shape == (4,) assert np_result.dtype == np.float64 - assert_array_equal(np_result, np.array([0., 1., 2., 3.])) + assert_array_equal(np_result, + np.array([0., 1., 2., 3., 0, 0, 0, 0, 0, 0, 0, 0, 0])) assert result.width == 3 assert result.max_degree == 2 @@ -76,9 +75,10 @@ def test_create_list_ints_no_ctype(): result = FreeTensor([0, 1, 2, 3], width=3, depth=2) np_result = np.array(result) - assert np_result.shape == (4,) + assert np_result.flags.owndata assert np_result.dtype == np.float64 - assert_array_equal(np_result, np.array([0., 1., 2., 3.])) + assert_array_equal(np_result, + np.array([0., 1., 2., 3., 0, 0, 0, 0, 0, 0, 0, 0, 0])) assert result.width == 3 assert result.max_degree == 2 @@ -111,7 +111,8 @@ def test_create_dlpack(): assert result.width == 3 assert result.max_degree == 2 - assert_array_equal(result, np.array([0., 1., 2., 3.])) + assert_array_equal(result, + np.array([0., 1., 2., 3., 0, 0, 0, 0, 0, 0, 0, 0, 0])) def test_create_buffer_doubles(): @@ -121,7 +122,8 @@ def test_create_buffer_doubles(): assert result.width == 3 assert result.max_degree == 2 - assert_array_equal(result, np.array([0., 1., 2., 3.])) + assert_array_equal(result, + np.array([0., 1., 2., 3., 0, 0, 0, 0, 0, 0, 0, 0, 0])) def test_create_buffer_floats(): @@ -131,7 +133,9 @@ def test_create_buffer_floats(): assert result.width == 3 assert result.max_degree == 2 - assert_array_equal(result, np.array([0., 1., 2., 3.], dtype=np.float32)) + assert_array_equal(result, + np.array([0., 1., 2., 3., 0, 0, 0, 0, 0, 0, 0, 0, 0], + dtype=np.float32)) def test_create_intv_pair(): @@ -181,7 +185,7 @@ def test_create_list_intv_pairs_dense(): assert result.width == 2 assert result.max_degree == 2 assert result.storage_type == roughpy.VectorType.DenseVector - assert_array_equal(result, np.array([0., 1., 2.])) + assert_array_equal(result, np.array([0., 1., 2., 0, 0, 0, 0])) def test_create_dict_args(): @@ -236,7 +240,7 @@ def test_create_FreeTensor_specified_width_incomplete_degree_range(rng, width): def test_FreeTensor_array_roundtrip(width, rdata, rtensor): assert rtensor.width == width - assert_array_equal(rdata, np.array(rtensor)) + assert_array_equal(rdata, np.array(rtensor)[:rdata.shape[0]]) def test_FreeTensor_repr(rng, width, depth, tensor_size): @@ -355,7 +359,7 @@ def test_FreeTensor_mul_single_letter(width): expected = FreeTensor(np.array( [1.0, 1.0] + [0.0] * (width - 1) + [0.5] + [0.0] * (width ** 2 - 1) + [ 1.0 / 6.0] + [0.0] * (width ** 3 - 1)), - width=width, depth=depth) + width=width, depth=depth) assert t.exp() == expected @@ -440,8 +444,8 @@ def test_antipode(width, depth, data1, coeff_type, vec_type): def test_free_tensor_pickle_roundtrip(): - - t = FreeTensor(np.array([1.0, 2.0, 3.0, 4.0]), width=3, depth=2, dtype=DPReal) + t = FreeTensor(np.array([1.0, 2.0, 3.0, 4.0]), width=3, depth=2, + dtype=DPReal) t2 = pickle.loads(pickle.dumps(t)) @@ -467,3 +471,57 @@ def test_ft_ctor_type_deduction(data, typ): f = FreeTensor(data, width=TYPE_DEDUCTION_WIDTH, depth=TYPE_DEDUCTION_DEPTH) assert f.dtype == typ + + +TO_ARRAY_WIDTH = 3 +TO_ARRAY_DEPTH = 3 + + +def to_array_tensor(): + ctx = roughpy.get_context(width=TO_ARRAY_WIDTH, depth=TO_ARRAY_DEPTH, + coeffs=roughpy.DPReal) + + yield FreeTensor(ctx=ctx) + yield FreeTensor([0.0], ctx=ctx) + + for depth in range(TO_ARRAY_WIDTH+1): + yield FreeTensor(np.zeros(ctx.tensor_size(depth)), ctx=ctx) + + # sparse + yield FreeTensor({0: 1.0, 1: 1.0}, ctx=ctx) + + +@pytest.mark.parametrize("tensor", to_array_tensor()) +def test_to_array(tensor): + ctx = roughpy.get_context(width=TO_ARRAY_WIDTH, depth=TO_ARRAY_DEPTH, + coeffs=roughpy.DPReal) + + array_tensor = np.array(tensor) + assert array_tensor.flags.owndata + assert_array_equal(array_tensor, + np.zeros(ctx.tensor_size(TO_ARRAY_WIDTH))) + + + + +def test_to_array_rational_coeffs(): + + ctx = roughpy.get_context(width=TO_ARRAY_WIDTH, depth=TO_ARRAY_DEPTH, coeffs=roughpy.Rational) + + tensor = FreeTensor([1]*(1+TO_ARRAY_WIDTH), ctx=ctx) + + array_tensor = np.array(tensor) + + assert array_tensor.dtype == object + +def test_to_array_rational_poly_coeffs(): + + ctx = roughpy.get_context(width=TO_ARRAY_WIDTH, depth=TO_ARRAY_DEPTH, coeffs=roughpy.RationalPoly) + + tensor = FreeTensor([1]*(1+TO_ARRAY_WIDTH), ctx=ctx) + + array_tensor = np.array(tensor) + + assert array_tensor.dtype == object + +