Skip to content

Commit

Permalink
Move sample_random into a method accessible outside
Browse files Browse the repository at this point in the history
  • Loading branch information
@ianhbell
ianhbell committed Jan 9, 2025
1 parent 1a6dfc8 commit fa55ead
Show file tree
Hide file tree
Showing 3 changed files with 151 additions and 40 deletions.
86 changes: 81 additions & 5 deletions doc/source/polygons/polygonvalidity.ipynb
View file
Original file line number Diff line number Diff line change
Expand Up @@ -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",
Expand All @@ -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"
]
Expand All @@ -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",
Expand All @@ -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": {
Expand Down
103 changes: 68 additions & 35 deletions include/teqpflsh/flasher.hpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -22,6 +22,71 @@ namespace teqpflsh {

using ArrayType = Eigen::ArrayXd;


/**
Derived from: https://stackoverflow.com/a/47418580
*/
template<typename Geo, typename Generator>
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 <class T, class = void>
struct has_operatorbracket1 : std::false_type {};

template <class T>
struct has_operatorbracket1 < T,
std::void_t<decltype(std::declval<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<typename Container>
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<double> 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<Container>::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
Expand Down Expand Up @@ -172,39 +237,7 @@ class QuadRegion2D{
/// area, and then randomly sampling a point inside a triangle
template<typename Container>
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<double> 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<Container>::value){
destx[rngcounter] = x;
desty[rngcounter] = y;
}
else{
destx(rngcounter) = x;
desty(rngcounter) = y;
}
}
return teqpflsh::sample_random<Container>(m_bounding_polygon.get(), Nsamples, destx, desty);
}

/**
Expand Down Expand Up @@ -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;
Expand Down
2 changes: 2 additions & 0 deletions interface/nanobind_interface.cpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -248,6 +248,8 @@ NB_MODULE(_teqpflsh_impl, m) {
.def("sample_gridded_w_tree", &QuadRegion2D::sample_gridded_w_tree<tensor1d>)
;

m.def("sample_random", &::sample_random<tensor1d>);

using env = QuadRegion2D::Envelope;
nb::class_<env>(m, "Envelope")
.def_ro("x_min", &env::x_min)
Expand Down

0 comments on commit fa55ead

Please sign in to comment.