diff --git a/doc/source/polygons/polygonvalidity.ipynb b/doc/source/polygons/polygonvalidity.ipynb index 547ad76..b6c0b51 100644 --- a/doc/source/polygons/polygonvalidity.ipynb +++ b/doc/source/polygons/polygonvalidity.ipynb @@ -14,19 +14,32 @@ "cell_type": "code", "execution_count": null, "id": "6b45a70c", - "metadata": {}, + "metadata": { + "editable": true, + "slideshow": { + "slide_type": "" + }, + "tags": [] + }, "outputs": [], "source": [ "import numpy as np\n", "import teqpflsh\n", - "import matplotlib.pyplot as plt" + "import matplotlib.pyplot as plt\n", + "teqpflsh.__version__" ] }, { "cell_type": "code", "execution_count": null, "id": "b6182393", - "metadata": {}, + "metadata": { + "editable": true, + "slideshow": { + "slide_type": "" + }, + "tags": [] + }, "outputs": [], "source": [ "ptr = teqpflsh.GeometryFactoryHolder()\n", @@ -42,8 +55,44 @@ }, { "cell_type": "markdown", - "id": "a8e9c1d5-3aa6-44b3-a766-9c413395e1c5", + "id": "ca4134dc-7c34-4894-a289-f51717086c14", "metadata": {}, + "source": [ + "The reason this matters in this context is that we need to be able to sample the domain randomly to generate guess values." + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "82c3f307-8322-4fa8-be80-f245ae866360", + "metadata": { + "editable": true, + "slideshow": { + "slide_type": "" + }, + "tags": [ + "raises-exception" + ] + }, + "outputs": [], + "source": [ + "# This doesn't work because the triangulation fails because the geometry is not valid\n", + "N = 10000\n", + "x = np.zeros((N, ))\n", + "y = np.zeros((N, ))\n", + "teqpflsh.sample_random(poly1, 10000, x, y)" + ] + }, + { + "cell_type": "markdown", + "id": "a8e9c1d5-3aa6-44b3-a766-9c413395e1c5", + "metadata": { + "editable": true, + "slideshow": { + "slide_type": "" + }, + "tags": [] + }, "source": [ "Now we need to break up the polygon into portions that are simple (non self-intersecting) with the MakeValid class of geos : https://libgeos.org/doxygen/classgeos_1_1operation_1_1valid_1_1MakeValid.html" ] @@ -52,7 +101,13 @@ "cell_type": "code", "execution_count": null, "id": "038ff472-6551-47c6-bb33-151f1ff4e483", - "metadata": {}, + "metadata": { + "editable": true, + "slideshow": { + "slide_type": "" + }, + "tags": [] + }, "outputs": [], "source": [ "simpl = poly1.make_valid()\n", @@ -66,6 +121,27 @@ " print(f'simple: {pI.isSimple}')\n", " print(f'valid: {pI.isValid}')" ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "30b47157-b3c5-4090-a4f4-10a642661cbb", + "metadata": { + "editable": true, + "slideshow": { + "slide_type": "" + }, + "tags": [] + }, + "outputs": [], + "source": [ + "# And now sampling works after forcing validity\n", + "N = 10000\n", + "x = np.zeros((N, ))\n", + "y = np.zeros((N, ))\n", + "teqpflsh.sample_random(simpl, 10000, x, y)\n", + "plt.plot(x, y, '.', ms=1);" + ] } ], "metadata": { diff --git a/include/teqpflsh/flasher.hpp b/include/teqpflsh/flasher.hpp index 340d5bc..c64e8ce 100644 --- a/include/teqpflsh/flasher.hpp +++ b/include/teqpflsh/flasher.hpp @@ -22,6 +22,71 @@ namespace teqpflsh { using ArrayType = Eigen::ArrayXd; + +/** + Derived from: https://stackoverflow.com/a/47418580 + */ +template +auto point_in_triangle(const Geo&tri, Generator& gen){ + std::uniform_real_distribution<> uni; + auto co = tri->getCoordinates(); + + auto r1 = uni(gen), r2 = uni(gen); + auto q = std::abs(r1 - r2); + auto s = q; + auto t = 0.5 * (r1 + r2 - q); + auto u = 1 - 0.5 * (q + r1 + r2); + auto x = s * co->getX(0) + t * co->getX(1) + u * co->getX(2); + auto y = s * co->getY(0) + t * co->getY(1) + u * co->getY(2); + return std::make_tuple(x, y); +} + +template +struct has_operatorbracket1 : std::false_type {}; + +template +struct has_operatorbracket1 < T, + std::void_t()[0])>> : + std::true_type {}; + +/// Sample a given number of points that are within the polygon by first doing triangulation of +/// the bounding polygon with triangulation, sampling the triangles according to their relative +/// area, and then randomly sampling a point inside a triangle +template +void sample_random(const Geometry* geo, std::size_t Nsamples, Container& destx, Container& desty){ + if(destx.size() != Nsamples){ throw std::invalid_argument("destx size does not equal Nsamples"); } + if(desty.size() != Nsamples){ throw std::invalid_argument("desty size does not equal Nsamples"); } + + auto tri = geos::triangulate::polygon::PolygonTriangulator::triangulate(geo); + + // Get the areas for each triangle + auto Ntri = tri->getNumGeometries(); + std::vector areas(Ntri); + for (auto igeo = 0U; igeo < Ntri; ++igeo){ + areas[igeo] = tri->getGeometryN(igeo)->getArea(); + } + // Prepare the weighting function + std::random_device rd; + std::mt19937 gen(rd()); + std::discrete_distribution<> d(areas.begin(), areas.end()); + + for (auto rngcounter = 0U; rngcounter < Nsamples; ++rngcounter){ + // Randomly select a triangle within the triangulation weighted according to + // its relative area. + auto selected_triangle = tri->getGeometryN(d(gen)); + // And select a point within this triangle + auto [x, y] = point_in_triangle(selected_triangle, gen); + if constexpr(has_operatorbracket1::value){ + destx[rngcounter] = x; + desty[rngcounter] = y; + } + else{ + destx(rngcounter) = x; + desty(rngcounter) = y; + } + } +} + /** This class is a generic type that represents a 2D region defined by its bounding envelope which is expected to be a closed & non-intersecting polygon. A quadtree can be constructed covering the @@ -172,39 +237,7 @@ class QuadRegion2D{ /// area, and then randomly sampling a point inside a triangle template void sample_random(std::size_t Nsamples, Container& destx, Container& desty){ - if(destx.size() != Nsamples){ throw std::invalid_argument("destx size does not equal Nsamples"); } - if(desty.size() != Nsamples){ throw std::invalid_argument("desty size does not equal Nsamples"); } - - - auto tri = do_fast_triangulation(); - - - // Get the areas for each triangle - auto Ntri = tri->getNumGeometries(); - std::vector areas(Ntri); - for (auto igeo = 0U; igeo < Ntri; ++igeo){ - areas[igeo] = tri->getGeometryN(igeo)->getArea(); - } - // Prepare the weighting function - std::random_device rd; - std::mt19937 gen(rd()); - std::discrete_distribution<> d(areas.begin(), areas.end()); - - for (auto rngcounter = 0U; rngcounter < Nsamples; ++rngcounter){ - // Randomly select a triangle within the triangulation weighted according to - // its relative area. - auto selected_triangle = tri->getGeometryN(d(gen)); - // And select a point within this triangle - auto [x, y] = point_in_triangle(selected_triangle, gen); - if constexpr(has_operatorbracket1::value){ - destx[rngcounter] = x; - desty[rngcounter] = y; - } - else{ - destx(rngcounter) = x; - desty(rngcounter) = y; - } - } + return teqpflsh::sample_random(m_bounding_polygon.get(), Nsamples, destx, desty); } /** @@ -783,11 +816,11 @@ class MainFlasher{ auto& r = optsoln.value(); // runtime_check_that(r.phases.size() > 0); T(i) = r.T_K; - rho(i) = r.rhobulk_molm3; +// rho(i) = r.rhobulk_molm3; q(i) = r.phases.front().qmolar; } }catch(std::exception&e){ - std::cout << e.what() << std::endl; +// std::cout << e.what() << std::endl; T(i) = -1; rho(i) = -1; q(i) = -1; diff --git a/interface/nanobind_interface.cpp b/interface/nanobind_interface.cpp index ab71438..f1a970e 100644 --- a/interface/nanobind_interface.cpp +++ b/interface/nanobind_interface.cpp @@ -248,6 +248,8 @@ NB_MODULE(_teqpflsh_impl, m) { .def("sample_gridded_w_tree", &QuadRegion2D::sample_gridded_w_tree) ; + m.def("sample_random", &::sample_random); + using env = QuadRegion2D::Envelope; nb::class_(m, "Envelope") .def_ro("x_min", &env::x_min)