// Copyright 2021 The Manifold Authors. // // Licensed under the Apache License, Version 2.0 (the "License"); // you may not use this file except in compliance with the License. // You may obtain a copy of the License at // // http://www.apache.org/licenses/LICENSE-2.0 // // Unless required by applicable law or agreed to in writing, software // distributed under the License is distributed on an "AS IS" BASIS, // WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied. // See the License for the specific language governing permissions and // limitations under the License. #pragma once #include "./parallel.h" #include "./utils.h" #include "./vec.h" #include "manifold/common.h" #include "manifold/optional_assert.h" namespace { template inline bool FirstFinite(T v) { return std::isfinite(v[0]); } template <> inline bool FirstFinite(double v) { return std::isfinite(v); } } // namespace namespace manifold { /** @ingroup Private */ class SparseIndices { // sparse indices where {p1: q1, p2: q2, ...} are laid out as // p1 q1 p2 q2 or q1 p1 q2 p2, depending on endianness // such that the indices are sorted by (p << 32) | q public: #if defined(__BYTE_ORDER) && __BYTE_ORDER == __BIG_ENDIAN || \ defined(__BIG_ENDIAN__) || defined(__ARMEB__) || defined(__THUMBEB__) || \ defined(__AARCH64EB__) || defined(_MIBSEB) || defined(__MIBSEB) || \ defined(__MIBSEB__) static constexpr size_t pOffset = 0; #elif defined(__BYTE_ORDER) && __BYTE_ORDER == __LITTLE_ENDIAN || \ defined(__LITTLE_ENDIAN__) || defined(__ARMEL__) || \ defined(__THUMBEL__) || defined(__AARCH64EL__) || defined(_MIPSEL) || \ defined(__MIPSEL) || defined(__MIPSEL__) || defined(__EMSCRIPTEN__) || \ defined(_WIN32) static constexpr size_t pOffset = 1; #else #error "unknown architecture" #endif static constexpr int64_t EncodePQ(int p, int q) { return (int64_t(p) << 32) | q; } SparseIndices() = default; SparseIndices(size_t size) { data_ = Vec(size * sizeof(int64_t)); } void Clear() { data_.clear(false); } void FromIndices(const std::vector& indices) { std::vector sizes; size_t total_size = 0; for (const auto& ind : indices) { sizes.push_back(total_size); total_size += ind.data_.size(); } data_ = Vec(total_size); for_each_n(ExecutionPolicy::Par, countAt(0), indices.size(), [&](size_t i) { std::copy(indices[i].data_.begin(), indices[i].data_.end(), data_.begin() + sizes[i]); }); } size_t size() const { return data_.size() / sizeof(int64_t); } Vec Copy(bool use_q) const { Vec out(size()); size_t offset = pOffset; if (use_q) offset = 1 - offset; const int* p = ptr(); for_each(autoPolicy(out.size()), countAt(0_uz), countAt(out.size()), [&](size_t i) { out[i] = p[i * 2 + offset]; }); return out; } void Sort() { VecView view = AsVec64(); stable_sort(view.begin(), view.end()); } void Resize(size_t size) { data_.resize(size * sizeof(int64_t), -1); } inline int& Get(size_t i, bool use_q) { if (use_q) return ptr()[2 * i + 1 - pOffset]; else return ptr()[2 * i + pOffset]; } inline int Get(size_t i, bool use_q) const { if (use_q) return ptr()[2 * i + 1 - pOffset]; else return ptr()[2 * i + pOffset]; } inline int64_t GetPQ(size_t i) const { VecView view = AsVec64(); return view[i]; } inline void Set(size_t i, int p, int q) { VecView view = AsVec64(); view[i] = EncodePQ(p, q); } inline void SetPQ(size_t i, int64_t pq) { VecView view = AsVec64(); view[i] = pq; } VecView AsVec64() { return VecView(reinterpret_cast(data_.data()), data_.size() / sizeof(int64_t)); } VecView AsVec64() const { return VecView( reinterpret_cast(data_.data()), data_.size() / sizeof(int64_t)); } VecView AsVec32() { return VecView(reinterpret_cast(data_.data()), data_.size() / sizeof(int32_t)); } VecView AsVec32() const { return VecView( reinterpret_cast(data_.data()), data_.size() / sizeof(int32_t)); } inline void Add(int p, int q, bool seq = false) { data_.extend(sizeof(int64_t), seq); Set(size() - 1, p, q); } void Unique() { Sort(); VecView view = AsVec64(); size_t newSize = unique(view.begin(), view.end()) - view.begin(); Resize(newSize); } size_t RemoveZeros(Vec& S) { DEBUG_ASSERT(S.size() == size(), userErr, "Different number of values than indicies!"); Vec new2Old(S.size()); sequence(new2Old.begin(), new2Old.end()); size_t size = copy_if(countAt(0_uz), countAt(S.size()), new2Old.begin(), [&S](const size_t i) { return S[i] != 0; }) - new2Old.begin(); new2Old.resize(size); Permute(S, new2Old); Vec tmp(std::move(data_)); Resize(size); gather(new2Old.begin(), new2Old.end(), reinterpret_cast(tmp.data()), reinterpret_cast(data_.data())); return size; } template size_t KeepFinite(Vec& v, Vec& x) { DEBUG_ASSERT(x.size() == size(), userErr, "Different number of values than indicies!"); Vec new2Old(v.size()); size_t size = copy_if(countAt(0_uz), countAt(v.size()), new2Old.begin(), [&v](size_t i) { return FirstFinite(v[i]); }) - new2Old.begin(); new2Old.resize(size); Permute(v, new2Old); Permute(x, new2Old); Vec tmp(std::move(data_)); Resize(size); gather(new2Old.begin(), new2Old.end(), reinterpret_cast(tmp.data()), reinterpret_cast(data_.data())); return size; } #ifdef MANIFOLD_DEBUG void Dump() const { std::cout << "SparseIndices = " << std::endl; const int* p = ptr(); for (size_t i = 0; i < size(); ++i) { std::cout << i << ", p = " << Get(i, false) << ", q = " << Get(i, true) << std::endl; } std::cout << std::endl; } #endif private: Vec data_; inline int* ptr() { return reinterpret_cast(data_.data()); } inline const int* ptr() const { return reinterpret_cast(data_.data()); } }; } // namespace manifold