From e56afd5412092e375133cb6014fd46d6770f78a0 Mon Sep 17 00:00:00 2001 From: Wenzel Jakob Date: Fri, 14 Feb 2020 14:08:17 +0100 Subject: [PATCH] Wrapped up distr_2d.h refactoring and fixed testcases --- .gitmodules | 3 - ext/nanogui | 2 +- ext/zeromq | 1 - include/mitsuba/core/distr_1d.h | 4 +- include/mitsuba/core/distr_2d.h | 606 +++++++++++++++-------------- setup.cfg | 5 +- src/bsdfs/CMakeLists.txt | 28 +- src/bsdfs/tests/test_diffuse.py | 11 +- src/conftest.py | 78 ++++ src/libcore/python/distr_2d_v.cpp | 2 +- src/libcore/tests/test_distr_1d.py | 57 ++- src/libcore/tests/test_distr_2d.py | 184 +++++++-- src/python/python/chi2.py | 21 +- src/python/python/test/__init__.py | 82 ---- 14 files changed, 605 insertions(+), 479 deletions(-) delete mode 160000 ext/zeromq create mode 100644 src/conftest.py diff --git a/.gitmodules b/.gitmodules index c4c0642d1..0e2d9d171 100644 --- a/.gitmodules +++ b/.gitmodules @@ -19,9 +19,6 @@ [submodule "ext/pugixml"] path = ext/pugixml url = https://github.com/mitsuba-renderer/pugixml -[submodule "ext/zeromq"] - path = ext/zeromq - url = https://github.com/mitsuba-renderer/zeromq [submodule "ext/asmjit"] path = ext/asmjit url = https://github.com/mitsuba-renderer/asmjit diff --git a/ext/nanogui b/ext/nanogui index 6b82fb18e..c8d66bba6 160000 --- a/ext/nanogui +++ b/ext/nanogui @@ -1 +1 @@ -Subproject commit 6b82fb18ed599017028b042ebb5aa831d0f93454 +Subproject commit c8d66bba6ed9244c38f34437dfe6c99f49e708db diff --git a/ext/zeromq b/ext/zeromq deleted file mode 160000 index d7c3a64ec..000000000 --- a/ext/zeromq +++ /dev/null @@ -1 +0,0 @@ -Subproject commit d7c3a64ecdda89209a49aee12c92c3a99da98c18 diff --git a/include/mitsuba/core/distr_1d.h b/include/mitsuba/core/distr_1d.h index 378d4ca51..dc52ab3ec 100644 --- a/include/mitsuba/core/distr_1d.h +++ b/include/mitsuba/core/distr_1d.h @@ -449,7 +449,7 @@ template struct ContinuousDistribution { Float y0 = gather(m_pdf, index, active), y1 = gather(m_pdf, index + 1u, active), - c0 = gather(m_cdf, index- 1u, active && index > 0); + c0 = gather(m_cdf, index - 1u, active && index > 0); value = (value - c0) * m_inv_interval_size; @@ -487,7 +487,7 @@ template struct ContinuousDistribution { Float y0 = gather(m_pdf, index, active), y1 = gather(m_pdf, index + 1u, active), - c0 = gather(m_cdf, index- 1u, active && index > 0); + c0 = gather(m_cdf, index - 1u, active && index > 0); value = (value - c0) * m_inv_interval_size; diff --git a/include/mitsuba/core/distr_2d.h b/include/mitsuba/core/distr_2d.h index 3e42ba0ee..f19ddb54e 100644 --- a/include/mitsuba/core/distr_2d.h +++ b/include/mitsuba/core/distr_2d.h @@ -65,13 +65,13 @@ template class Distribution2D { Throw("Distribution2D(): input array resolution must be >= 2!"); // The linear interpolant has 'size-1' patches - ScalarVector2u n_patches = size - 1u; + ScalarVector2u n_patches = size - 1; m_patch_size = 1.f / n_patches; m_inv_patch_size = n_patches; // Dependence on additional parameters - m_slices = 1u; + m_slices = 1; for (int i = (int) Dimension - 1; i >= 0; --i) { if (param_res[i] < 1) Throw("Distribution2D(): parameter resolution must be >= 1!"); @@ -101,7 +101,8 @@ template class Distribution2D { UInt32 param_index = math::find_interval( (uint32_t) m_param_values[dim].size(), [&](UInt32 idx) ENOKI_INLINE_LAMBDA { - return gather(m_param_values[dim], idx, active) < param[dim]; + return gather(m_param_values[dim], idx, active) < + param[dim]; }); Float p0 = gather(m_param_values[dim], param_index, active), @@ -118,7 +119,7 @@ template class Distribution2D { ENOKI_MARK_USED(param); ENOKI_MARK_USED(param_weight); ENOKI_MARK_USED(active); - return 0.f; + return 0u; } } @@ -205,7 +206,7 @@ class Hierarchical2D : public Distribution2D { * still be sampled (proportionally), but returned density values * will reflect the unnormalized values. * - * If \c build_hierarchy is set to \c false, the implementation will not + * If \c enable_sampling is set to \c false, the implementation will not * construct the hierarchy needed for sample warping, which saves memory * in case this functionality is not needed (e.g. if only the interpolation * in \c eval() is used). In this case, \c sample() and \c invert() @@ -217,18 +218,18 @@ class Hierarchical2D : public Distribution2D { const std::array ¶m_res = { }, const std::array ¶m_values = { }, bool normalize = true, - bool build_hierarchy = true) + bool enable_sampling = true) : Base(size, param_res, param_values) { // The linear interpolant has 'size-1' patches - ScalarVector2u n_patches = size - 1u; + ScalarVector2u n_patches = size - 1; // Keep track of the dependence on additional parameters (optional) uint32_t max_level = math::log2i_ceil(hmax(n_patches)); - m_max_patch_index = n_patches - 1u; + m_max_patch_index = n_patches - 1; - if (!build_hierarchy) { + if (!enable_sampling) { m_levels.reserve(1); m_levels.emplace_back(size, m_slices); @@ -335,20 +336,20 @@ class Hierarchical2D : public Distribution2D { // Fetch values from next MIP level UInt32 offset_i = level.index(offset) + slice_offset * level.size; - Float v00 = level.template lookup<>(offset_i, m_param_strides, - param_weight, active); + Float v00 = level.lookup(offset_i, m_param_strides, + param_weight, active); offset_i += 1u; - Float v10 = level.template lookup<>(offset_i, m_param_strides, - param_weight, active); + Float v10 = level.lookup(offset_i, m_param_strides, + param_weight, active); offset_i += 1u; - Float v01 = level.template lookup<>(offset_i, m_param_strides, - param_weight, active); + Float v01 = level.lookup(offset_i, m_param_strides, + param_weight, active); offset_i += 1u; - Float v11 = level.template lookup<>(offset_i, m_param_strides, - param_weight, active); + Float v11 = level.lookup(offset_i, m_param_strides, + param_weight, active); // Avoid issues with roundoff error sample = clamp(sample, 0.f, 1.f); @@ -378,19 +379,17 @@ class Hierarchical2D : public Distribution2D { offset.x() + offset.y() * level0.width + slice_offset * level0.size; // Fetch corners of bilinear patch - Float v00 = level0.template lookup<>(offset_i, m_param_strides, - param_weight, active); + Float v00 = level0.lookup(offset_i, m_param_strides, + param_weight, active); - Float v10 = level0.template lookup<>(offset_i + 1, m_param_strides, - param_weight, active); + Float v10 = level0.lookup(offset_i + 1, m_param_strides, + param_weight, active); - Float v01 = level0.template lookup<>(offset_i + level0.width, - m_param_strides, param_weight, - active); + Float v01 = level0.lookup(offset_i + level0.width, m_param_strides, + param_weight, active); - Float v11 = level0.template lookup<>(offset_i + level0.width + 1, - m_param_strides, param_weight, - active); + Float v11 = level0.lookup(offset_i + level0.width + 1, m_param_strides, + param_weight, active); Float pdf; std::tie(sample, pdf) = @@ -403,7 +402,8 @@ class Hierarchical2D : public Distribution2D { } /// Inverse of the mapping implemented in \c sample() - std::pair invert(Point2f sample, const Float *param = nullptr, + std::pair invert(Point2f sample, + const Float *param = nullptr, Mask active = true) const { MTS_MASK_ARGUMENT(active); @@ -424,19 +424,17 @@ class Hierarchical2D : public Distribution2D { UInt32 offset_i = offset.x() + offset.y() * level0.width + slice_offset * level0.size; - Float v00 = level0.template lookup<>(offset_i, m_param_strides, - param_weight, active); + Float v00 = level0.lookup(offset_i, m_param_strides, + param_weight, active); - Float v10 = level0.template lookup<>(offset_i + 1, m_param_strides, - param_weight, active); + Float v10 = level0.lookup(offset_i + 1, m_param_strides, + param_weight, active); - Float v01 = level0.template lookup<>(offset_i + level0.width, - m_param_strides, param_weight, - active); + Float v01 = level0.lookup(offset_i + level0.width, m_param_strides, + param_weight, active); - Float v11 = level0.template lookup<>(offset_i + level0.width + 1, - m_param_strides, param_weight, - active); + Float v11 = level0.lookup(offset_i + level0.width + 1, m_param_strides, + param_weight, active); sample -= Point2f(Point2i(offset)); @@ -450,20 +448,20 @@ class Hierarchical2D : public Distribution2D { // Fetch values from next MIP level offset_i = level.index(offset & ~1u) + slice_offset * level.size; - v00 = level.template lookup<>(offset_i, m_param_strides, - param_weight, active); + v00 = level.lookup(offset_i, m_param_strides, + param_weight, active); offset_i += 1u; - v10 = level.template lookup<>(offset_i, m_param_strides, - param_weight, active); + v10 = level.lookup(offset_i, m_param_strides, + param_weight, active); offset_i += 1u; - v01 = level.template lookup<>(offset_i, m_param_strides, - param_weight, active); + v01 = level.lookup(offset_i, m_param_strides, + param_weight, active); offset_i += 1u; - v11 = level.template lookup<>(offset_i, m_param_strides, - param_weight, active); + v11 = level.lookup(offset_i, m_param_strides, + param_weight, active); Mask x_mask = neq(offset.x() & 1u, 0u), y_mask = neq(offset.y() & 1u, 0u); @@ -512,26 +510,24 @@ class Hierarchical2D : public Distribution2D { UInt32 offset_i = offset.x() + offset.y() * level0.width + slice_offset * level0.size; - Float v00 = level0.template lookup<>(offset_i, m_param_strides, - param_weight, active); + Float v00 = level0.lookup(offset_i, m_param_strides, + param_weight, active); - Float v10 = level0.template lookup<>(offset_i + 1, m_param_strides, - param_weight, active); + Float v10 = level0.lookup(offset_i + 1, m_param_strides, + param_weight, active); - Float v01 = level0.template lookup<>(offset_i + level0.width, - m_param_strides, param_weight, - active); + Float v01 = level0.lookup(offset_i + level0.width, m_param_strides, + param_weight, active); - Float v11 = level0.template lookup<>(offset_i + level0.width + 1, - m_param_strides, param_weight, - active); + Float v11 = level0.lookup(offset_i + level0.width + 1, m_param_strides, + param_weight, active); return warp::square_to_bilinear_pdf(v00, v10, v01, v11, pos); } std::string to_string() const { std::ostringstream oss; - oss << "Hierarchical2D<" << Dimension << ">[" << std::endl + oss << "Hierarchical2D" << Dimension << "[" << std::endl << " size = [" << m_levels[0].width << ", " << m_levels[0].size / m_levels[0].width << "]," << std::endl << " levels = " << m_levels.size() << "," << std::endl; @@ -614,11 +610,7 @@ class Hierarchical2D : public Distribution2D { } else { ENOKI_MARK_USED(param_strides); ENOKI_MARK_USED(param_weight); - - if constexpr (is_scalar_v) - return gather(data, i0); - else - return gather(data, i0, active); + return gather(data, i0, active); } } }; @@ -655,8 +647,8 @@ class Hierarchical2D : public Distribution2D { * evaluating the distribution for in-between parameter values. * * \remark The Python API exposes explicitly instantiated versions of this - * class named \c Marginal2DDiscrete0 to \c Marginal2DDiscrete3 and - * \c Marginal2DContinuous0 to \c Marginal2DContinuous3 and for data that depends + * class named \c MarginalDiscrete2D0 to \c MarginalDiscrete2D3 and \c + * MarginalContinuous2D0 to \c MarginalContinuous2D3 for data that depends * on 0 to 3 parameters. */ template @@ -688,7 +680,7 @@ class Marginal2D : public Distribution2D { * still be sampled (proportionally), but returned density values * will reflect the unnormalized values. * - * If \c build_cdf is set to \c false, the implementation will not + * If \c enable_sampling is set to \c false, the implementation will not * construct the cdf needed for sample warping, which saves memory in case * this functionality is not needed (e.g. if only the interpolation in \c * eval() is used). @@ -697,94 +689,95 @@ class Marginal2D : public Distribution2D { const ScalarVector2u &size, const std::array ¶m_res = { }, const std::array ¶m_values = { }, - bool normalize = true, bool build_cdf = true) + bool normalize = true, bool enable_sampling = true) : Base(size, param_res, param_values), m_size(size), m_normalized(normalize) { - uint32_t w = m_size.x(), - h = m_size.y(), - n_values = w * h, - n_marginal = h - 1, - n_conditional; + uint32_t w = m_size.x(), + h = m_size.y(), + n_data = w * h, + n_marg = h - 1, + n_cond = (w - 1) * (Continuous ? h : (h - 1)); - if constexpr (Continuous) - n_conditional = h * (w - 1); - else - n_conditional = (h - 1) * (w - 1); + double scale_x = .5 / (w - 1), + scale_y = .5 / (h - 1); - m_data = empty(m_slices * n_values); + m_data = empty(m_slices * n_data); m_data.managed(); - if (build_cdf) { - m_marg_cdf = empty(m_slices * n_marginal); + if (enable_sampling) { + m_marg_cdf = empty(m_slices * n_marg); m_marg_cdf.managed(); - m_cond_cdf = empty(m_slices * n_conditional); + m_cond_cdf = empty(m_slices * n_cond); m_cond_cdf.managed(); - ScalarFloat *marginal_cdf = m_marg_cdf.data(), - *conditional_cdf = m_cond_cdf.data(), + ScalarFloat *marg_cdf = m_marg_cdf.data(), + *cond_cdf = m_cond_cdf.data(), *data_out = m_data.data(); - double scale = hprod(1.0 / (m_size - 1)) - * (Continuous ? .5 : 0.25); + std::unique_ptr cond_cdf_sum(new double[h]); for (uint32_t slice = 0; slice < m_slices; ++slice) { - /* The continuous / discrete cases require different - tables with marginal/conditional probabilities */ - double sum = 0.0; + ScalarFloat norm = 1.f; + + /* The marginal/probability distribution computation + differs for the Continuous=false/true cases */ if constexpr (Continuous) { // Construct conditional CDF for (uint32_t y = 0; y < h; ++y) { + double accum = 0.0; uint32_t i = y * w, j = y * (w - 1); for (uint32_t x = 0; x < w - 1; ++x, ++i, ++j) { - sum += scale * ((double) data[i] + - (double) data[i + 1]); - conditional_cdf[j] = (ScalarFloat) sum; + accum += scale_x * ((double) data[i] + + (double) data[i + 1]); + cond_cdf[j] = (ScalarFloat) accum; } + cond_cdf_sum[y] = accum; } // Construct marginal CDF - sum = 0.0; + double accum = 0.0; for (uint32_t y = 0; y < h - 1; ++y) { - sum += .5 * ((double) conditional_cdf[(y + 1) * (w - 1) - 1] + - (double) conditional_cdf[(y + 2) * (w - 1) - 1]); - marginal_cdf[y] = (ScalarFloat) sum; + accum += scale_y * (cond_cdf_sum[y] + cond_cdf_sum[y + 1]); + marg_cdf[y] = (ScalarFloat) accum; } + + if (normalize) + norm = ScalarFloat(1.0 / accum); } else { + double scale = scale_x * scale_y; + // Construct conditional CDF for (uint32_t y = 0; y < h - 1; ++y) { + double accum = 0.0; uint32_t i = y * w, j = y * (w - 1); for (uint32_t x = 0; x < w - 1; ++x, ++i, ++j) { - sum += scale * ((double) data[i] + - (double) data[i + 1] + - (double) data[i + w] + - (double) data[i + w + 1]); - conditional_cdf[j] = (ScalarFloat) sum; + accum += scale * ((double) data[i] + + (double) data[i + 1] + + (double) data[i + w] + + (double) data[i + w + 1]); + cond_cdf[j] = (ScalarFloat) accum; } + cond_cdf_sum[y] = accum; } // Construct marginal CDF - sum = 0.0; + double accum = 0.0; for (uint32_t y = 0; y < h - 1; ++y) { - sum += (double) conditional_cdf[(y + 1) * (w - 1) - 1]; - marginal_cdf[y] = (ScalarFloat) sum; + accum += cond_cdf_sum[y]; + marg_cdf[y] = (ScalarFloat) accum; } + + if (normalize) + norm = ScalarFloat(1.0 / accum); } - // Normalize CDFs and PDF (if requested) - ScalarFloat norm = 1.f; - if (normalize) - norm = ScalarFloat(1.0 / sum); - - for (size_t i = 0; i < n_conditional; ++i) - *conditional_cdf++ *= norm; - for (size_t i = 0; i < n_marginal; ++i) - *marginal_cdf++ *= norm; - for (size_t i = 0; i < n_values; ++i) + for (size_t i = 0; i < n_cond; ++i) + *cond_cdf++ *= norm; + for (size_t i = 0; i < n_marg; ++i) + *marg_cdf++ *= norm; + for (size_t i = 0; i < n_data; ++i) *data_out++ = *data++ * norm; - - marginal_cdf += m_size.y(); - conditional_cdf += n_values; } } else { ScalarFloat *data_out = m_data.data(); @@ -797,18 +790,16 @@ class Marginal2D : public Distribution2D { for (uint32_t y = 0; y < h - 1; ++y) { size_t i = y * w; for (uint32_t x = 0; x < w - 1; ++x, ++i) { - ScalarFloat v00 = data[i], - v10 = data[i + 1], - v01 = data[i + w], - v11 = data[i + w + 1], - avg = .25f * (v00 + v10 + v01 + v11); - sum += (double) avg; + sum += (double) data[i] + + (double) data[i + 1] + + (double) data[i + w] + + (double) data[i + w + 1]; } } - norm = ScalarFloat(1.0 / sum); + norm = ScalarFloat(1.0 / (scale_x * scale_y * sum)); } - for (uint32_t k = 0; k < n_values; ++k) + for (uint32_t k = 0; k < n_data; ++k) *data_out++ = *data++ * norm; } } @@ -823,6 +814,7 @@ class Marginal2D : public Distribution2D { std::pair sample(const Point2f &sample, const Float *param = nullptr, Mask active = true) const { + Assert(!m_marg_cdf.empty(), "Marginal2D::sample(): enable_sampling=false!"); if constexpr (Continuous) return sample_continuous(sample, param, active); @@ -834,6 +826,8 @@ class Marginal2D : public Distribution2D { std::pair invert(const Point2f &sample, const Float *param = nullptr, Mask active = true) const { + Assert(!m_marg_cdf.empty(), "Marginal2D::invert(): enable_sampling=false!"); + if constexpr (Continuous) return invert_continuous(sample, param, active); else @@ -865,21 +859,21 @@ class Marginal2D : public Distribution2D { if (Dimension != 0) index += slice_offset * size; - Float v00 = lookup<>(m_data.data(), index, - size, param_weight, active), - v10 = lookup<>(m_data.data() + 1, index, - size, param_weight, active), - v01 = lookup<>(m_data.data() + m_size.x(), index, - size, param_weight, active), - v11 = lookup<>(m_data.data() + m_size.x() + 1, index, - size, param_weight, active); + Float v00 = lookup(m_data.data(), index, + size, param_weight, active), + v10 = lookup(m_data.data() + 1, index, + size, param_weight, active), + v01 = lookup(m_data.data() + m_size.x(), index, + size, param_weight, active), + v11 = lookup(m_data.data() + m_size.x() + 1, index, + size, param_weight, active); return warp::square_to_bilinear_pdf(v00, v10, v01, v11, pos); } std::string to_string() const { std::ostringstream oss; - oss << "Marginal2D<" << Dimension << ">[" << std::endl + oss << "Marginal2D" << Dimension << "[" << std::endl << " size = " << m_size << "," << std::endl; if (Dimension > 0) { oss << " param_size = ["; @@ -935,9 +929,9 @@ class Marginal2D : public Distribution2D { MTS_MASK_ARGUMENT(active); // Size of a slice of various tables (conditional/marginal/data) - uint32_t size_cond = hprod(m_size - 1), - size_marg = m_size.y() - 1, - size_data = hprod(m_size); + uint32_t n_cond = hprod(m_size - 1), + n_marg = m_size.y() - 1, + n_data = hprod(m_size); /// Find offset and interpolation weights wrt. conditional parameters Float param_weight[2 * DimensionInt]; @@ -947,69 +941,69 @@ class Marginal2D : public Distribution2D { sample = clamp(sample, math::Epsilon, math::OneMinusEpsilon); /// Multiply by last entry of marginal CDF if the data is not normalized - UInt32 offset_marg = slice_offset * size_marg; + UInt32 offset_marg = slice_offset * n_marg; + + auto fetch_marginal = [&](UInt32 idx, Mask mask) + ENOKI_INLINE_LAMBDA -> Float { + return lookup(m_marg_cdf.data(), offset_marg + idx, + n_marg, param_weight, mask); + }; + if (!m_normalized) - sample.y() *= lookup<>(m_marg_cdf.data(), - offset_marg + size_marg - 1, - size_marg, param_weight, active); + sample.y() *= fetch_marginal(n_marg - 1, active); // Sample the row from the marginal distribution UInt32 row = enoki::binary_search( - 0u, size_marg - 1u, [&](UInt32 idx) ENOKI_INLINE_LAMBDA { - return lookup<>(m_marg_cdf.data(), offset_marg + idx, - size_marg, param_weight, - active) < sample.y(); + 0u, n_marg - 1, [&](UInt32 idx) ENOKI_INLINE_LAMBDA { + return fetch_marginal(idx, active) < sample.y(); }); // Re-scale uniform variate - Float row_cdf_0 = lookup<>(m_marg_cdf.data(), offset_marg + row - 1, - size_marg, param_weight, active && row > 0), - row_cdf_1 = lookup<>(m_marg_cdf.data(), offset_marg + row, - size_marg, param_weight, active); + Float row_cdf_0 = fetch_marginal(row - 1, active && row > 0), + row_cdf_1 = fetch_marginal(row, active); sample.y() -= row_cdf_0; masked(sample.y(), neq(row_cdf_1, row_cdf_0)) /= row_cdf_1 - row_cdf_0; /// Multiply by last entry of conditional CDF - UInt32 offset_cond = - slice_offset * size_cond + row * (m_size.x() - 1); - sample.x() *= lookup<>(m_cond_cdf.data(), - offset_cond + m_size.x() - 2, - size_cond, param_weight, active); + UInt32 offset_cond = slice_offset * n_cond + row * (m_size.x() - 1); + sample.x() *= lookup(m_cond_cdf.data(), + offset_cond + m_size.x() - 2, + n_cond, param_weight, active); // Sample the column from the conditional distribution UInt32 col = enoki::binary_search( - 0u, size_cond - 1u, [&](UInt32 idx) ENOKI_INLINE_LAMBDA { - return lookup<>(m_cond_cdf.data(), - offset_cond + idx, size_cond, - param_weight, active) < sample.x(); + 0u, m_size.x() - 2, [&](UInt32 idx) ENOKI_INLINE_LAMBDA { + return lookup(m_cond_cdf.data(), + offset_cond + idx, n_cond, + param_weight, active) < sample.x(); }); // Re-scale uniform variate - Float col_cdf_0 = lookup<>(m_cond_cdf.data(), offset_cond + col - 1, - size_cond, param_weight, active && col > 0), - col_cdf_1 = lookup<>(m_cond_cdf.data(), offset_cond + col, - size_cond, param_weight, active); + Float col_cdf_0 = lookup(m_cond_cdf.data(), offset_cond + col - 1, + n_cond, param_weight, active && col > 0), + col_cdf_1 = lookup(m_cond_cdf.data(), offset_cond + col, + n_cond, param_weight, active); sample.x() -= col_cdf_0; masked(sample.x(), neq(col_cdf_1, col_cdf_0)) /= col_cdf_1 - col_cdf_0; // Sample a position on the bilinear patch - UInt32 offset_data = slice_offset * size_data + row * m_size.x() + col; - Float v00 = lookup<>(m_data.data(), offset_data, size_data, - param_weight, active), - v10 = lookup<>(m_data.data() + 1, offset_data, size_data, - param_weight, active), - v01 = lookup<>(m_data.data() + m_size.x(), offset_data, - size_data, param_weight, active), - v11 = lookup<>(m_data.data() + m_size.x() + 1, offset_data, - size_data, param_weight, active); + UInt32 offset_data = slice_offset * n_data + row * m_size.x() + col; + Float v00 = lookup(m_data.data(), offset_data, n_data, + param_weight, active), + v10 = lookup(m_data.data() + 1, offset_data, n_data, + param_weight, active), + v01 = lookup(m_data.data() + m_size.x(), offset_data, + n_data, param_weight, active), + v11 = lookup(m_data.data() + m_size.x() + 1, offset_data, + n_data, param_weight, active); Float pdf; std::tie(sample, pdf) = warp::square_to_bilinear(v00, v10, v01, v11, sample); - return { (Point2u(col, row) + sample) * m_patch_size, pdf }; + return { (Point2i(Point2u(col, row)) + sample) * m_patch_size, pdf }; } MTS_INLINE @@ -1019,9 +1013,9 @@ class Marginal2D : public Distribution2D { MTS_MASK_ARGUMENT(active); // Size of a slice of various tables (conditional/marginal/data) - uint32_t size_cond = hprod(m_size - 1), - size_marg = m_size.y() - 1, - size_data = hprod(m_size); + uint32_t n_cond = hprod(m_size - 1), + n_marg = m_size.y() - 1, + n_data = hprod(m_size); /// Find offset and interpolation weights wrt. conditional parameters Float param_weight[2 * DimensionInt]; @@ -1033,44 +1027,43 @@ class Marginal2D : public Distribution2D { // Fetch values at corners of bilinear patch sample *= m_inv_patch_size; Point2u offset = min(Point2u(Point2i(sample)), m_size - 2u); - UInt32 index = offset.x() + offset.y() * m_size.x() + slice_offset * size_data; + UInt32 index = offset.x() + offset.y() * m_size.x() + slice_offset * n_data; sample -= Point2f(Point2i(offset)); - Float v00 = lookup<>(m_data.data(), index, - size_data, param_weight, active), - v10 = lookup<>(m_data.data() + 1, index, - size_data, param_weight, active), - v01 = lookup<>(m_data.data() + m_size.x(), index, - size_data, param_weight, active), - v11 = lookup<>(m_data.data() + m_size.x() + 1, index, - size_data, param_weight, active); + Float v00 = lookup(m_data.data(), index, + n_data, param_weight, active), + v10 = lookup(m_data.data() + 1, index, + n_data, param_weight, active), + v01 = lookup(m_data.data() + m_size.x(), index, + n_data, param_weight, active), + v11 = lookup(m_data.data() + m_size.x() + 1, index, + n_data, param_weight, active); Float pdf; std::tie(sample, pdf) = warp::bilinear_to_square(v00, v10, v01, v11, sample); - UInt32 cond_offset = slice_offset * size_cond + offset.y() * (m_size.x() - 1), - marg_offset = slice_offset * size_marg; + UInt32 offset_cond = slice_offset * n_cond + offset.y() * (m_size.x() - 1), + offset_marg = slice_offset * n_marg; - Float row_cdf_0 = lookup<>(m_marg_cdf.data(), marg_offset + offset.y() - 1, - size_marg, param_weight, active && offset.y() > 0), - row_cdf_1 = lookup<>(m_marg_cdf.data(), marg_offset + offset.y(), - size_marg, param_weight, active), - col_cdf_0 = lookup<>(m_cond_cdf.data(), cond_offset + offset.x() - 1, - size_cond, param_weight, active && offset.x() > 0), - col_cdf_1 = lookup<>(m_cond_cdf.data(), cond_offset + offset.x(), - size_cond, param_weight, active); + Float row_cdf_0 = lookup(m_marg_cdf.data(), offset_marg + offset.y() - 1, + n_marg, param_weight, active && offset.y() > 0), + row_cdf_1 = lookup(m_marg_cdf.data(), offset_marg + offset.y(), + n_marg, param_weight, active), + col_cdf_0 = lookup(m_cond_cdf.data(), offset_cond + offset.x() - 1, + n_cond, param_weight, active && offset.x() > 0), + col_cdf_1 = lookup(m_cond_cdf.data(), offset_cond + offset.x(), + n_cond, param_weight, active); - sample.x() = fmadd(sample.x(), col_cdf_1 - col_cdf_0, col_cdf_0); - sample.y() = fmadd(sample.y(), row_cdf_1 - row_cdf_0, row_cdf_0); + sample.x() = lerp(col_cdf_0, col_cdf_1, sample.x()); + sample.y() = lerp(row_cdf_0, row_cdf_1, sample.y()); - sample.x() /= lookup<>(m_cond_cdf.data(), - cond_offset + m_size.x() - 2, - size_cond, param_weight, active); + sample.x() /= lookup(m_cond_cdf.data(), + offset_cond + m_size.x() - 2, + n_cond, param_weight, active); if (!m_normalized) - sample.y() /= lookup<>(m_marg_cdf.data(), - marg_offset + m_size.y() - 2, - size_marg, param_weight, active); + sample.y() /= lookup(m_marg_cdf.data(), offset_marg + n_marg - 1, + n_marg, param_weight, active); return { sample, pdf }; } @@ -1081,90 +1074,92 @@ class Marginal2D : public Distribution2D { Mask active) const { MTS_MASK_ARGUMENT(active); + // Size of a slice of various tables (conditional/marginal/data) + uint32_t n_cond = m_size.y() * (m_size.x() - 1), + n_marg = m_size.y() - 1, + n_data = hprod(m_size); + + /// Find offset and interpolation weights wrt. conditional parameters Float param_weight[2 * DimensionInt]; UInt32 slice_offset = interpolate_weights(param, param_weight, active); ENOKI_MARK_USED(slice_offset); // Avoid degeneracies on the domain boundary - sample = clamp(sample, math::Epsilon, math::OneMinusEpsilon); + sample = clamp(sample, math::Epsilon, + math::OneMinusEpsilon); // Sample the row first - UInt32 marg_offset = 0; - if constexpr (Dimension != 0) - marg_offset = slice_offset * m_size.y(); + UInt32 offset_marg = slice_offset * n_marg; - auto fetch_marginal = [&](UInt32 idx) ENOKI_INLINE_LAMBDA -> Float { - return lookup<>(m_marg_cdf.data(), marg_offset + idx, - m_size.y(), param_weight, active); + auto fetch_marginal = [&](UInt32 idx, Mask mask) + ENOKI_INLINE_LAMBDA -> Float { + return lookup(m_marg_cdf.data(), offset_marg + idx, + n_marg, param_weight, mask); }; - UInt32 row = math::find_interval( - m_size.y(), + if (!m_normalized) + sample.y() *= fetch_marginal(n_marg - 1, active); + + UInt32 row = enoki::binary_search( + 0u, n_marg - 1, [&](UInt32 idx) ENOKI_INLINE_LAMBDA { - return fetch_marginal(idx) < sample.y(); + return fetch_marginal(idx, active) < sample.y(); }); - sample.y() -= fetch_marginal(row); - - uint32_t slice_size = hprod(m_size); - UInt32 cond_offset = row * m_size.x(); - if constexpr (Dimension != 0) - cond_offset += slice_offset * slice_size; + /// Subtract the marginal CDF value up to the current interval + sample.y() -= fetch_marginal(row - 1, active && row > 0); - Float r0 = lookup<>(m_cond_cdf.data(), - cond_offset + m_size.x() - 1, slice_size, - param_weight, active), - r1 = lookup<>(m_cond_cdf.data(), - cond_offset + (m_size.x() * 2 - 1), slice_size, - param_weight, active); + UInt32 offset_cond = slice_offset * n_cond + row * (m_size.x() - 1); - Mask is_const = abs(r0 - r1) < 1e-4f * (r0 + r1); - sample.y() = select(is_const, 2.f * sample.y(), - r0 - safe_sqrt(r0 * r0 - 2.f * sample.y() * (r0 - r1))); - sample.y() /= select(is_const, r0 + r1, r0 - r1); + /// Look up conditional CDF values of surrounding rows for x == 1 + Float r0 = lookup(m_cond_cdf.data() + 1 * (m_size.x() - 1) - 1, + offset_cond, n_cond, param_weight, active), + r1 = lookup(m_cond_cdf.data() + 2 * (m_size.x() - 1) - 1, + offset_cond, n_cond, param_weight, active); - // Sample the column next - sample.x() *= (1.f - sample.y()) * r0 + sample.y() * r1; + sample.y() = sample_segment(sample.y(), m_inv_patch_size.y(), r0, r1); - auto fetch_conditional = [&](UInt32 idx) ENOKI_INLINE_LAMBDA -> Float { - Float v0 = lookup<>(m_cond_cdf.data(), cond_offset + idx, - slice_size, param_weight, active), - v1 = lookup<>(m_cond_cdf.data() + m_size.x(), - cond_offset + idx, slice_size, param_weight, active); + // Multiply sample.x() by the integrated density along the 'x' axis + sample.x() *= lerp(r0, r1, sample.y()); - return (1.f - sample.y()) * v0 + sample.y() * v1; + // Sample the column next + auto fetch_conditional = [&](UInt32 idx, Mask mask) + ENOKI_INLINE_LAMBDA -> Float { + idx += offset_cond; + Float v0 = lookup(m_cond_cdf.data(), + idx, n_cond, param_weight, mask), + v1 = lookup(m_cond_cdf.data() + m_size.x() - 1, + idx, n_cond, param_weight, mask); + return lerp(v0, v1, sample.y()); }; - UInt32 col = math::find_interval( - m_size.x(), + UInt32 col = enoki::binary_search( + 0, m_size.x() - 1, [&](UInt32 idx) ENOKI_INLINE_LAMBDA { - return fetch_conditional(idx) < sample.x(); + return fetch_conditional(idx, active) < sample.x(); }); - sample.x() -= fetch_conditional(col); + /// Subtract the CDF value up to the current interval + sample.x() -= fetch_conditional(col - 1, active && col > 0); - cond_offset += col; + UInt32 offset_data = slice_offset * n_data + row * m_size.x() + col; - Float v00 = lookup<>(m_data.data(), cond_offset, slice_size, - param_weight, active), - v10 = lookup<>(m_data.data() + 1, cond_offset, slice_size, - param_weight, active), - v01 = lookup<>(m_data.data() + m_size.x(), cond_offset, - slice_size, param_weight, active), - v11 = lookup<>(m_data.data() + m_size.x() + 1, cond_offset, - slice_size, param_weight, active), - c0 = fmadd((1.f - sample.y()), v00, sample.y() * v01), - c1 = fmadd((1.f - sample.y()), v10, sample.y() * v11); + Float v00 = lookup(m_data.data(), offset_data, n_data, + param_weight, active), + v10 = lookup(m_data.data() + 1, offset_data, n_data, + param_weight, active), + v01 = lookup(m_data.data() + m_size.x(), offset_data, + n_data, param_weight, active), + v11 = lookup(m_data.data() + m_size.x() + 1, offset_data, + n_data, param_weight, active), + c0 = lerp(v00, v01, sample.y()), + c1 = lerp(v10, v11, sample.y()); - is_const = abs(c0 - c1) < 1e-4f * (c0 + c1); - sample.x() = select(is_const, 2.f * sample.x(), - c0 - safe_sqrt(c0 * c0 - 2.f * sample.x() * (c0 - c1))); - Float divisor = select(is_const, c0 + c1, c0 - c1); - masked(sample.x(), neq(divisor, 0.f)) /= divisor; + sample.x() = sample_segment(sample.x(), m_inv_patch_size.x(), c0, c1); return { - (Point2u(col, row) + sample) * m_patch_size, - fmadd(1.f - sample.x(), c0, sample.x() * c1) * hprod(m_inv_patch_size) + (Point2i(Point2u(col, row)) + sample) * m_patch_size, + lerp(c0, c1, sample.x()) }; } @@ -1174,6 +1169,12 @@ class Marginal2D : public Distribution2D { Mask active) const { MTS_MASK_ARGUMENT(active); + // Size of a slice of various tables (conditional/marginal/data) + uint32_t n_cond = m_size.y() * (m_size.x() - 1), + n_marg = m_size.y() - 1, + n_data = hprod(m_size); + + /// Find offset and interpolation weights wrt. conditional parameters Float param_weight[2 * DimensionInt]; UInt32 slice_offset = interpolate_weights(param, param_weight, active); @@ -1185,60 +1186,75 @@ class Marginal2D : public Distribution2D { Point2u pos = min(Point2u(Point2i(sample)), m_size - 2u); sample -= Point2f(Point2i(pos)); - UInt32 offset = pos.x() + pos.y() * m_size.x(); - uint32_t slice_size = hprod(m_size); - if constexpr (Dimension != 0) - offset += slice_offset * slice_size; + UInt32 offset_data = + slice_offset * n_data + pos.y() * m_size.x() + pos.x(); // Invert the X component - Float v00 = lookup<>(m_data.data(), offset, slice_size, - param_weight, active), - v10 = lookup<>(m_data.data() + 1, offset, slice_size, - param_weight, active), - v01 = lookup<>(m_data.data() + m_size.x(), offset, slice_size, - param_weight, active), - v11 = lookup<>(m_data.data() + m_size.x() + 1, offset, slice_size, - param_weight, active); + Float v00 = lookup(m_data.data(), offset_data, n_data, + param_weight, active), + v10 = lookup(m_data.data() + 1, offset_data, n_data, + param_weight, active), + v01 = lookup(m_data.data() + m_size.x(), offset_data, + n_data, param_weight, active), + v11 = lookup(m_data.data() + m_size.x() + 1, offset_data, + n_data, param_weight, active); - Point2f w1 = sample, w0 = 1.f - w1; + Float c0 = lerp(v00, v01, sample.y()), + c1 = lerp(v10, v11, sample.y()), + pdf = lerp(c0, c1, sample.x()); - Float c0 = fmadd(w0.y(), v00, w1.y() * v01), - c1 = fmadd(w0.y(), v10, w1.y() * v11), - pdf = fmadd(w0.x(), c0, w1.x() * c1); + sample.x() = invert_segment(sample.x(), m_patch_size.x(), c0, c1); - sample.x() *= c0 + .5f * sample.x() * (c1 - c0); + UInt32 offset_cond = + slice_offset * n_cond + pos.y() * (m_size.x() - 1); + + auto fetch_conditional = [&](UInt32 idx, Mask mask) + ENOKI_INLINE_LAMBDA -> Float { + idx += offset_cond; + Float v0 = lookup(m_cond_cdf.data(), + idx, n_cond, param_weight, mask), + v1 = lookup(m_cond_cdf.data() + m_size.x() - 1, + idx, n_cond, param_weight, mask); + return lerp(v0, v1, sample.y()); + }; - Float v0 = lookup<>(m_cond_cdf.data(), offset, - slice_size, param_weight, active), - v1 = lookup<>(m_cond_cdf.data() + m_size.x(), - offset, slice_size, param_weight, active); + sample.x() += fetch_conditional(pos.x() - 1, active && pos.x() > 0); - sample.x() += (1.f - sample.y()) * v0 + sample.y() * v1; + Float r0 = lookup(m_cond_cdf.data() + 1 * (m_size.x() - 1) - 1, + offset_cond, n_cond, param_weight, active), + r1 = lookup(m_cond_cdf.data() + 2 * (m_size.x() - 1) - 1, + offset_cond, n_cond, param_weight, active); - offset = pos.y() * m_size.x(); - if (Dimension != 0) - offset += slice_offset * slice_size; + sample.x() /= lerp(r0, r1, sample.y()); - Float r0 = lookup<>(m_cond_cdf.data(), - offset + m_size.x() - 1, slice_size, - param_weight, active), - r1 = lookup<>(m_cond_cdf.data(), - offset + (m_size.x() * 2 - 1), slice_size, - param_weight, active); + // Invert the Y component + sample.y() = invert_segment(sample.y(), m_patch_size.y(), r0, r1); - sample.x() /= (1.f - sample.y()) * r0 + sample.y() * r1; + UInt32 offset_marg = slice_offset * n_marg; + sample.y() += lookup(m_marg_cdf.data(), offset_marg + pos.y() - 1, + m_size.y(), param_weight, active && pos.y() > 0); - // Invert the Y component - sample.y() *= r0 + .5f * sample.y() * (r1 - r0); + if (!m_normalized) + sample.y() /= lookup(m_marg_cdf.data(), offset_marg + n_marg - 1, + n_marg, param_weight, active); - offset = pos.y(); - if (Dimension != 0) - offset += slice_offset * m_size.y(); + return { sample, pdf }; + } - sample.y() += lookup<>(m_marg_cdf.data(), offset, - m_size.y(), param_weight, active); + MTS_INLINE Float sample_segment(Float sample, ScalarFloat inv_width, + Float v0, Float v1) const { + Mask non_const = abs(v0 - v1) > 1e-4f * (v0 + v1); + Float divisor = select(non_const, v0 - v1, v0 + v1); + sample *= 2.f * inv_width; + masked(sample, non_const) = + v0 - safe_sqrt(sqr(v0) + sample * (v1 - v0)); + masked(sample, neq(divisor, 0.f)) /= divisor; + return sample; + } - return { sample, pdf * hprod(m_inv_patch_size) }; + MTS_INLINE Float invert_segment(Float sample, ScalarFloat width, + Float v0, Float v1) const { + return sample * lerp(v0, v1, .5f * sample) * width; } protected: /// Resolution of the discretized density function diff --git a/setup.cfg b/setup.cfg index 16d8f7ffa..210bcdab6 100644 --- a/setup.cfg +++ b/setup.cfg @@ -1,8 +1,9 @@ [tool:pytest] minversion = 3.0 -norecursedirs = ext ext_build build build-debug CMakeFiles dist include +norecursedirs = ext ext_build build build-debug CMakeFiles dist include .git python_paths = dist dist/python + [pycodestyle] -# E402 module level import not at top of file +# E402 module level import not at top of file (needed for Mitsuba variants) ignore = E402 diff --git a/src/bsdfs/CMakeLists.txt b/src/bsdfs/CMakeLists.txt index cdc43c5d5..b25a0aea5 100644 --- a/src/bsdfs/CMakeLists.txt +++ b/src/bsdfs/CMakeLists.txt @@ -1,20 +1,20 @@ set(MTS_PLUGIN_PREFIX "bsdfs") -add_plugin(blendbsdf blendbsdf.cpp) -add_plugin(conductor conductor.cpp) -add_plugin(dielectric dielectric.cpp) -add_plugin(diffuse diffuse.cpp) -add_plugin(mask mask.cpp) -add_plugin(measured measured.cpp) -add_plugin(null null.cpp) -add_plugin(plastic plastic.cpp) -add_plugin(roughconductor roughconductor.cpp) +add_plugin(blendbsdf blendbsdf.cpp) +add_plugin(conductor conductor.cpp) +add_plugin(dielectric dielectric.cpp) +add_plugin(diffuse diffuse.cpp) +add_plugin(mask mask.cpp) +add_plugin(measured measured.cpp) +add_plugin(null null.cpp) +add_plugin(plastic plastic.cpp) +add_plugin(roughconductor roughconductor.cpp) add_plugin(roughdielectric roughdielectric.cpp) -add_plugin(roughplastic roughplastic.cpp) -add_plugin(thindielectric thindielectric.cpp) -add_plugin(twosided twosided.cpp) -add_plugin(polarizer polarizer.cpp) -add_plugin(retarder retarder.cpp) +add_plugin(roughplastic roughplastic.cpp) +add_plugin(thindielectric thindielectric.cpp) +add_plugin(twosided twosided.cpp) +add_plugin(polarizer polarizer.cpp) +add_plugin(retarder retarder.cpp) # Register the test directory add_tests(${CMAKE_CURRENT_SOURCE_DIR}/tests) diff --git a/src/bsdfs/tests/test_diffuse.py b/src/bsdfs/tests/test_diffuse.py index 4b9a32dea..3e73ef047 100644 --- a/src/bsdfs/tests/test_diffuse.py +++ b/src/bsdfs/tests/test_diffuse.py @@ -43,17 +43,16 @@ def test02_eval_pdf(variant_scalar): def test03_chi2(variant_packet): from mitsuba.python.chi2 import BSDFAdapter, ChiSquareTest, SphericalDomain - from mitsuba.core import ScalarBoundingBox2f sample_func, pdf_func = BSDFAdapter("diffuse", '') chi2 = ChiSquareTest( - domain = SphericalDomain(), - sample_func = sample_func, - pdf_func = pdf_func, - sample_dim = 3 + domain=SphericalDomain(), + sample_func=sample_func, + pdf_func=pdf_func, + sample_dim=3 ) result = chi2.run(0.1) # chi2._dump_tables() - assert result \ No newline at end of file + assert result diff --git a/src/conftest.py b/src/conftest.py new file mode 100644 index 000000000..fc4a4db77 --- /dev/null +++ b/src/conftest.py @@ -0,0 +1,78 @@ +# This file rewrites exceptions caught by PyTest and makes traces involving +# pybind11 classes more readable. It e.g. replaces "" by "allclose" which is arguably just +# as informative and much more compact. + +import pytest +import re + +re1 = re.compile(r'') +re2 = re.compile(r']*>') + + +def patch_line(s): + s = re1.sub(r'\1', s) + s = re2.sub(r'\1', s) + return s + + +def patch_tb(tb): # tb: ReprTraceback + for entry in tb.reprentries: + entry.lines = [patch_line(l) for l in entry.lines] + + +@pytest.hookimpl(hookwrapper=True) +def pytest_runtest_makereport(item, call): + outcome = yield + rep = outcome.get_result() + if rep.outcome == 'failed': + for entry in rep.longrepr.chain: + patch_tb(entry[0]) + return rep + + +def generate_fixture(variant): + @pytest.fixture(scope='module') + def fixture(): + try: + import mitsuba + mitsuba.set_variant(variant) + except Exception: + pytest.skip('Mitsuba variant "%s" is not enabled!' % variant) + globals()['variant_' + variant] = fixture + + +for variant in ['scalar_rgb', 'scalar_spectral', + 'scalar_mono_polarized', 'packet_rgb']: + generate_fixture(variant) +del generate_fixture + + +@pytest.fixture(params=['packet_rgb', 'gpu_rgb']) +def variants_vec_rgb(request): + try: + import mitsuba + mitsuba.set_variant(request.param) + except Exception: + pytest.skip('Mitsuba variant "%s" is not enabled!' % request.param) + return request.param + + +@pytest.fixture(params=['scalar_rgb', 'packet_rgb']) +def variants_cpu_rgb(request): + try: + import mitsuba + mitsuba.set_variant(request.param) + except Exception: + pytest.skip('Mitsuba variant "%s" is not enabled!' % request.param) + return request.param + + +@pytest.fixture(params=['scalar_rgb', 'packet_rgb', 'gpu_rgb']) +def variants_all_rgb(request): + try: + import mitsuba + mitsuba.set_variant(request.param) + except Exception: + pytest.skip('Mitsuba variant "%s" is not enabled!' % request.param) + return request.param diff --git a/src/libcore/python/distr_2d_v.cpp b/src/libcore/python/distr_2d_v.cpp index 3829b2512..4d7d75e4a 100644 --- a/src/libcore/python/distr_2d_v.cpp +++ b/src/libcore/python/distr_2d_v.cpp @@ -49,7 +49,7 @@ template auto bind_warp(py::module &m, if constexpr (Warp::Dimension == 0) warp.def(std::move(constructor), "data"_a, "param_values"_a = py::list(), "normalize"_a = true, - "build_hierarchy"_a = true, doc_constructor); + "enable_sampling"_a = true, doc_constructor); else warp.def(std::move(constructor), "data"_a, "param_values"_a, "normalize"_a = true, "build_hierarchy"_a = true, diff --git a/src/libcore/tests/test_distr_1d.py b/src/libcore/tests/test_distr_1d.py index d32fa4e7a..5adb29cd4 100644 --- a/src/libcore/tests/test_distr_1d.py +++ b/src/libcore/tests/test_distr_1d.py @@ -1,11 +1,9 @@ -import mitsuba import pytest import enoki as ek -from mitsuba.python.test import variant_packet - -def test01_discr_empty(variant_packet): +def test01_discr_empty(variant_packet_rgb): + # Test that operations involving the empty distribution throw from mitsuba.core import DiscreteDistribution d = DiscreteDistribution() @@ -16,7 +14,8 @@ def test01_discr_empty(variant_packet): assert 'empty distribution' in str(excinfo.value) -def test02_discr_zero_prob(variant_packet): +def test02_discr_zero_prob(variant_packet_rgb): + # Test that operations involving zero probability mass throw from mitsuba.core import DiscreteDistribution with pytest.raises(RuntimeError) as excinfo: @@ -24,7 +23,8 @@ def test02_discr_zero_prob(variant_packet): assert "no probability mass found" in str(excinfo.value) -def test03_discr_neg_prob(variant_packet): +def test03_discr_neg_prob(variant_packet_rgb): + # Test that operations involving negative probability mass throw from mitsuba.core import DiscreteDistribution with pytest.raises(RuntimeError) as excinfo: @@ -32,7 +32,8 @@ def test03_discr_neg_prob(variant_packet): assert "entries must be non-negative" in str(excinfo.value) -def test04_discr_basic(variant_packet): +def test04_discr_basic(variant_packet_rgb): + # Validate discrete distribution cdf/pmf against hand-computed reference from mitsuba.core import DiscreteDistribution, Float x = DiscreteDistribution([1, 3, 2]) @@ -63,7 +64,8 @@ def test04_discr_basic(variant_packet): assert ek.allclose(x.normalization(), 1.0 / 3.0) -def test05_discr_sample(variant_packet): +def test05_discr_sample(variant_packet_rgb): + # Validate discrete distribution sampling against hand-computed reference from mitsuba.core import DiscreteDistribution, Float eps = 1e-7 @@ -100,7 +102,8 @@ def test05_discr_sample(variant_packet): ) -def test06_discr_bruteforce(variant_packet): +def test06_discr_bruteforce(variant_packet_rgb): + # Brute force validation of discrete distribution sampling from mitsuba.core import DiscreteDistribution, Float, PCG32, UInt64 rng = PCG32(initseq=UInt64.arange(50)) @@ -117,10 +120,12 @@ def test06_discr_bruteforce(variant_packet): z = ek.gather(ddistr.cdf(), y - 1, y > 0) x *= ddistr.sum() + # Did we sample the right interval? assert ek.all((x > z) | (ek.eq(x, 0) & (x >= z))) -def test07_discr_leading_trailing_zeros(variant_packet): +def test07_discr_leading_trailing_zeros(variant_packet_rgb): + # Check that sampling still works when there are zero-valued buckets from mitsuba.core import DiscreteDistribution x = DiscreteDistribution([0, 0, 1, 0, 1, 0, 0, 0]) index, pmf = x.sample_pmf([-100, 0, 0.5, 0.5 + 1e-6, 1, 100]) @@ -128,7 +133,8 @@ def test07_discr_leading_trailing_zeros(variant_packet): assert pmf == [.5] * 6 -def test08_cont_empty(variant_packet): +def test08_cont_empty(variant_packet_rgb): + # Test that operations involving the empty distribution throw from mitsuba.core import ContinuousDistribution d = ContinuousDistribution() @@ -144,7 +150,8 @@ def test08_cont_empty(variant_packet): assert 'needs at least two entries' in str(excinfo.value) -def test09_cont_empty_invalid_range(variant_packet): +def test09_cont_empty_invalid_range(variant_packet_rgb): + # Test that invalid range specifications throw an exception from mitsuba.core import ContinuousDistribution with pytest.raises(RuntimeError) as excinfo: @@ -156,7 +163,8 @@ def test09_cont_empty_invalid_range(variant_packet): assert 'invalid range' in str(excinfo.value) -def test10_cont_zero_prob(variant_packet): +def test10_cont_zero_prob(variant_packet_rgb): + # Test that operations involving zero probability mass throw from mitsuba.core import ContinuousDistribution with pytest.raises(RuntimeError) as excinfo: @@ -164,7 +172,8 @@ def test10_cont_zero_prob(variant_packet): assert "no probability mass found" in str(excinfo.value) -def test11_cont_neg_prob(variant_packet): +def test11_cont_neg_prob(variant_packet_rgb): + # Test that operations involving negative probability mass throw from mitsuba.core import ContinuousDistribution with pytest.raises(RuntimeError) as excinfo: @@ -172,7 +181,8 @@ def test11_cont_neg_prob(variant_packet): assert "entries must be non-negative" in str(excinfo.value) -def test12_cont_eval(variant_packet): +def test12_cont_eval(variant_packet_rgb): + # Test continuous 1D distribution pdf/cdf against hand-computed reference from mitsuba.core import ContinuousDistribution d = ContinuousDistribution([2, 3], [1, 2]) eps = 1e-6 @@ -198,7 +208,8 @@ def test12_cont_eval(variant_packet): ) -def test13_cont_func(variant_packet): +def test13_cont_func(variant_packet_rgb): + # Test continuous 1D distribution integral against analytic result from mitsuba.core import ContinuousDistribution, Float x = ek.linspace(Float, -2, 2, 513) @@ -210,7 +221,8 @@ def test13_cont_func(variant_packet): assert ek.allclose(d.sample([0, 0.5, 1]), [-2, 0, 2]) -def test14_irrcont_empty(variant_packet): +def test14_irrcont_empty(variant_packet_rgb): + # Test that operations involving the empty distribution throw from mitsuba.core import IrregularContinuousDistribution d = IrregularContinuousDistribution() @@ -230,7 +242,8 @@ def test14_irrcont_empty(variant_packet): assert 'size mismatch' in str(excinfo.value) -def test15_irrcont_empty_invalid_range(variant_packet): +def test15_irrcont_empty_invalid_range(variant_packet_rgb): + # Test that invalid range specifications throw an exception from mitsuba.core import IrregularContinuousDistribution with pytest.raises(RuntimeError) as excinfo: @@ -242,7 +255,8 @@ def test15_irrcont_empty_invalid_range(variant_packet): assert 'strictly increasing' in str(excinfo.value) -def test16_irrcont_zero_prob(variant_packet): +def test16_irrcont_zero_prob(variant_packet_rgb): + # Test that operations involving the empty distribution throw from mitsuba.core import IrregularContinuousDistribution with pytest.raises(RuntimeError) as excinfo: @@ -250,7 +264,8 @@ def test16_irrcont_zero_prob(variant_packet): assert "no probability mass found" in str(excinfo.value) -def test17_irrcont_neg_prob(variant_packet): +def test17_irrcont_neg_prob(variant_packet_rgb): + # Test that operations involving negative probability mass throw from mitsuba.core import IrregularContinuousDistribution with pytest.raises(RuntimeError) as excinfo: @@ -258,7 +273,7 @@ def test17_irrcont_neg_prob(variant_packet): assert "entries must be non-negative" in str(excinfo.value) -def test18_irrcont_simple_function(variant_packet): +def test18_irrcont_simple_function(variant_packet_rgb): # Reference from Mathematica from mitsuba.core import IrregularContinuousDistribution, Float diff --git a/src/libcore/tests/test_distr_2d.py b/src/libcore/tests/test_distr_2d.py index 8ac1ac848..6edce7bcb 100644 --- a/src/libcore/tests/test_distr_2d.py +++ b/src/libcore/tests/test_distr_2d.py @@ -3,19 +3,14 @@ import enoki as ek import numpy as np -from mitsuba.python.test import variant_scalar, variant_packet, variants_vec +@pytest.mark.parametrize("normalize", [True, False]) +@pytest.mark.parametrize("warp", ['Hierarchical2D0', 'MarginalDiscrete2D0']) +def test01_sample_inverse_discrete(variant_scalar_rgb, warp, normalize): + # Spot checks of Hierarchial2D0/MarginalDiscrete2D0 vs Mathematica -@pytest.fixture(params=['Hierarchical2D', 'MarginalDiscrete2D', - 'MarginalContinuous2D']) -def warps(request): - return [getattr(mitsuba.core, request.param + str(i)) for i in range(4)] - - -def test01_sample_inverse_simple(variant_scalar, warps): - # Spot checks of Hierarchical2D0, MarginalDiscrete2D0 - if 'Continuous' in warps[0].__name__: - return + def allclose(a, b): + return ek.allclose(a, b, atol=1e-6) # Mismatched X/Y resolution, odd number of columns ref = np.array([[1, 2, 5], @@ -24,45 +19,152 @@ def test01_sample_inverse_simple(variant_scalar, warps): # Normalized discrete PDF corresponding to the two patches intg = np.array([19, 16]) / 35 - distr = warps[0](ref) + warp = getattr(mitsuba.core, warp) + distr = warp(ref, normalize=normalize) + s = 35 / 8.0 if not normalize else 1 # Check if we can reach the corners and transition between patches - assert ek.allclose(distr.sample([0, 0]), [[0, 0], 8.0 / 35.0], atol=1e-6) - assert ek.allclose(distr.sample([1, 1]), [[1, 1], 16.0 / 35.0], atol=1e-6) - assert ek.allclose(distr.sample([intg[0], 0]), [[0.5, 0], 16.0 / 35.0], - atol=1e-6) - - assert ek.allclose(distr.invert([0, 0]), [[0, 0], 8.0 / 35.0], atol=1e-6) - assert ek.allclose(distr.invert([1, 1]), [[1, 1], 16.0 / 35.0], atol=1e-6) - assert ek.allclose(distr.invert([0.5, 0]), [[intg[0], 0], 16.0 / 35.0], - atol=1e-6) + assert allclose(distr.sample([0, 0]), [[0, 0], s * 8.0 / 35.0]) + assert allclose(distr.sample([1, 1]), [[1, 1], s * 16.0 / 35.0]) + assert allclose(distr.sample([intg[0], 0]), [[0.5, 0], s * 16.0 / 35.0]) + assert allclose(distr.invert([0, 0]), [[0, 0], s * 8.0 / 35.0]) + assert allclose(distr.invert([1, 1]), [[1, 1], s * 16.0 / 35.0]) + assert allclose(distr.invert([0.5, 0]), [[intg[0], 0], s * 16.0 / 35.0]) from mitsuba.core.warp import bilinear_to_square # Check if we can sample a specific position in each patch sample, pdf = bilinear_to_square(1, 2, 9, 7, [0.4, 0.3]) sample.x *= intg[0] - pdf *= 8.0 / 35.0 - assert ek.allclose(distr.sample(sample), [[0.2, 0.3], pdf]) - assert ek.allclose(distr.invert([0.2, 0.3]), [sample, pdf]) - assert ek.allclose(distr.eval([0.2, 0.3]), pdf) + pdf *= 8.0 / 35.0 * s + assert allclose(distr.sample(sample), [[0.2, 0.3], pdf]) + assert allclose(distr.invert([0.2, 0.3]), [sample, pdf]) + assert allclose(distr.eval([0.2, 0.3]), pdf) sample, pdf = bilinear_to_square(2, 5, 7, 2, [0.4, 0.3]) sample.x = sample.x * intg[1] + intg[0] - pdf *= 8.0 / 35.0 - assert ek.allclose(distr.sample(sample), [[0.7, 0.3], pdf]) - assert ek.allclose(distr.invert([0.7, 0.3]), [sample, pdf]) - assert ek.allclose(distr.eval([0.7, 0.3]), pdf+1) - -def test01_sample_simple(variant_scalar, warps): - if 'Continuous' in warps[0].__name__: - return - -# Test inverse<->sample, eval -# Chi2 test pow2, npow2 -# both hi res and low res -# Test interpolation with and without normalization -# Test sampling with and without interpolation -# Test combinations of multiple distributins at boundary and in the middle -# Does invert still work for non-normalized data? + pdf *= 8.0 / 35.0 * s + assert allclose(distr.sample(sample), [[0.7, 0.3], pdf]) + assert allclose(distr.invert([0.7, 0.3]), [sample, pdf]) + assert allclose(distr.eval([0.7, 0.3]), pdf) + +@pytest.mark.parametrize("normalize", [True, False]) +def test02_sample_inverse_continuous(variant_scalar_rgb, normalize): + # Spot checks of MarginalContinuous2D0 vs Mathematica + + from mitsuba.core import MarginalContinuous2D0 + + # Nonuniform case + ref = np.array([[1, 2, 5], + [9, 7, 2], + [1, 0, 5]]) + + distr = MarginalContinuous2D0(ref, normalize=normalize) + s = 1 if normalize else 4.125 + assert ek.allclose(distr.eval([0, 0]), s / 4.125) + ref = [[.162433, .330829], s * 1.44807] + assert ek.allclose(distr.sample([0.2, 0.3]), ref) + assert ek.allclose(distr.invert(ref[0]), [[0.2, 0.3], ref[1]]) + ref = [[.8277673, .8235800], s * 0.8326217] + assert ek.allclose(distr.sample([0.8, 0.9]), ref, atol=1e-5) + assert ek.allclose(distr.invert(ref[0]), [[0.8, 0.9], ref[1]]) + + # Uniform case + ref = np.array([[1, 1, 1], + [1, 1, 1], + [1, 1, 1]]) + + distr = MarginalContinuous2D0(ref) + assert ek.allclose(distr.eval([0, 0]), 1.0) + assert ek.allclose(distr.sample([0.2, 0.3]), + [[.2, .3], 1.0]) + assert ek.allclose(distr.invert([0.2, 0.3]), + [[.2, .3], 1.0]) + assert ek.allclose(distr.sample([0.8, 0.9]), + [[.8, .9], 1.0]) + assert ek.allclose(distr.invert([0.8, 0.9]), + [[.8, .9], 1.0]) + + +all_warps = [w + str(i) for w in + ['Hierarchical2D', 'MarginalDiscrete2D', 'MarginalContinuous2D'] + for i in range(4)] + + +@pytest.mark.parametrize("normalize", [True, False]) +@pytest.mark.parametrize("warp", all_warps) +def test03_forward_inverse_nd(variant_scalar_rgb, warp, normalize): + # Check that forward + inverse give the identity. Test for + # all supported variants of N-dimensional warps (N = 0..3) + + cls = getattr(mitsuba.core, warp) + ndim = int(warp[-1]) + 2 + np.random.seed(all_warps.index(warp)) + + for i in range(10): + shape = np.random.randint(2, 8, ndim) + param_res = [ + sorted(np.random.rand(s)) for s in shape[:-2] + ] + values = np.random.rand(*shape) * 10 + if i == 9: + values = np.ones(shape) + instance = cls(values, param_res, normalize=normalize) + + for j in range(10): + p_i = np.random.rand(2) + p_o, pdf = instance.sample(p_i) + assert ek.allclose(pdf, instance.eval(p_o), atol=1e-4) + p_i_2, pdf2 = instance.invert(p_o) + assert ek.allclose(pdf, pdf2, atol=1e-4) + assert ek.allclose(p_i_2, p_i, atol=1e-4) + + +@pytest.mark.parametrize("attempt", list(range(10))) +@pytest.mark.parametrize("warp", all_warps) +def test04_chi2(variant_packet_rgb, warp, attempt): + # Chi^2 test to make sure that the warping schemes are producing + # the right distribution. Test for all supported variants of + # N-dimensional warps (N = 0..3) + + from mitsuba.core import ScalarBoundingBox2f + from mitsuba.python.chi2 import ChiSquareTest, PlanarDomain + + cls = getattr(mitsuba.core, warp) + ndim = int(warp[-1]) + 2 + np.random.seed(all_warps.index(warp) * 10 + attempt) + + shape = np.random.randint(2, 8, ndim) + param_res = [ + sorted(np.random.rand(s)) for s in shape[:-2] + ] + + if attempt == 9: + values = np.ones(shape) + else: + values = np.random.rand(*shape) * 10 + + instance = cls(values, param_res) + + for j in range(10 if ndim > 2 else 1): + param = [ + ek.lerp( + param_res[i][0], + param_res[i][-1], + np.random.rand() + ) for i in range(0, ndim - 2) + ] + + chi2 = ChiSquareTest( + domain=PlanarDomain(ScalarBoundingBox2f(0, 1)), + sample_func=lambda p: instance.sample(p, param=param)[0], + pdf_func=lambda p: instance.eval(p, param=param), + sample_dim=2, + res=31, + sample_count=100000 + ) + + assert chi2.run( + test_count=11 * len(all_warps) + ) diff --git a/src/python/python/chi2.py b/src/python/python/chi2.py index f7460f62b..e2624dacc 100644 --- a/src/python/python/chi2.py +++ b/src/python/python/chi2.py @@ -213,7 +213,7 @@ def tabulate_pdf(self): self.pdf *= self.sample_count * ek.hprod(self.bounds.extents()) - def run(self, significance_level, test_count=1): + def run(self, significance_level=0.01, test_count=1, quiet=False): """ Run the Chi^2 test @@ -278,7 +278,7 @@ def run(self, significance_level, test_count=1): # Probability of observing a test statistic at least as # extreme as the one here assuming that the distributions match - self.p_value = 1-rlgamma(dof/2, chi2val/2) + self.p_value = 1 - rlgamma(dof / 2, chi2val / 2) # Apply the Šidák correction term, since we'll be conducting multiple # independent hypothesis tests. This accounts for the fact that the @@ -288,23 +288,26 @@ def run(self, significance_level, test_count=1): (1.0 - significance_level) ** (1.0 / test_count) if self.fail: - self._dump_tables() self._log('Not running the test for reasons listed above. Target ' 'density and histogram were written to "chi2_data.py') - return False + result = False elif self.p_value < significance_level \ or not ek.isfinite(self.p_value): - self._dump_tables() self._log('***** Rejected ***** the null hypothesis (p-value = %f,' ' significance level = %f). Target density and histogram' ' were written to "chi2_data.py".' % (self.p_value, significance_level)) - return False + result = False else: self._log('Accepted the null hypothesis (p-value = %f, ' 'significance level = %f)' % (self.p_value, significance_level)) - return True + result = True + if not quiet: + print(self.messages) + if not result: + self._dump_tables() + return result def _dump_tables(self): with open("chi2_data.py", "w") as f: @@ -548,6 +551,4 @@ def my_pdf(p): sample_dim=2 ) - chi2.run(0.1) - print(chi2.messages) - chi2._dump_tables() + chi2.run(0.01) diff --git a/src/python/python/test/__init__.py b/src/python/python/test/__init__.py index def24d7fe..e69de29bb 100644 --- a/src/python/python/test/__init__.py +++ b/src/python/python/test/__init__.py @@ -1,82 +0,0 @@ -import pytest -import mitsuba - -@pytest.fixture() -def variant_scalar(): - mitsuba.set_variant('scalar_rgb') - - -@pytest.fixture() -def variant_packet(): - try: - mitsuba.set_variant('packet_rgb') - except: - pytest.skip("packet_rgb mode not enabled") - - -@pytest.fixture() -def variant_spectral(): - try: - mitsuba.set_variant('scalar_spectral') - except: - pytest.skip("scalar_spectral mode not enabled") - - -@pytest.fixture() -def variant_polarized(): - try: - mitsuba.set_variant('scalar_spectral_polarized') - except: - pytest.skip("scalar_spectral_polarized mode not enabled") - - -@pytest.fixture() -def variant_mono_polarized(): - try: - mitsuba.set_variant('scalar_mono_polarized') - except: - pytest.skip("scalar_mono_polarized mode not enabled") - - -@pytest.fixture() -def variant_packet_spectral(): - try: - mitsuba.set_variant('packet_spectral') - except: - pytest.skip("packet_spectral mode not enabled") - - -@pytest.fixture(params = ['packet_rgb', 'gpu_rgb']) -def variants_vec(request): - try: - mitsuba.set_variant(request.param) - except: - pytest.skip("%s mode not enabled" % request.param) - return request.param - - -@pytest.fixture(params = ['scalar_rgb', 'packet_rgb']) -def variants_cpu_rgb(request): - try: - mitsuba.set_variant(request.param) - except: - pytest.skip("%s mode not enabled" % request.param) - return request.param - - -@pytest.fixture(params = ['scalar_rgb', 'packet_rgb', 'gpu_rgb']) -def variants_all_rgb(request): - try: - mitsuba.set_variant(request.param) - except: - pytest.skip("%s mode not enabled" % request.param) - return request.param - - -@pytest.fixture(params = mitsuba.variants()) -def variants_all(request): - try: - mitsuba.set_variant(request.param) - except: - pytest.skip("%s mode not enabled" % request.param) - return request.param \ No newline at end of file