diff --git a/.DS_Store b/.DS_Store new file mode 100644 index 000000000..b45760b19 Binary files /dev/null and b/.DS_Store differ diff --git a/.gitignore b/.gitignore index 7fcdd7723..cf3ef1dcc 100644 --- a/.gitignore +++ b/.gitignore @@ -13,10 +13,10 @@ site *.grib *.gif *.html +*.req +.mypy_cache example_eo example_mri -.mypy_cache -*.req polytope_python.egg-info new_test_venv_polytope new_venv_polytope @@ -26,3 +26,12 @@ new_updated_numpy_venv newest-polytope-venv serializedTree new_polytope_venv + +newest-polytope-venv +*.pstats +*.profile +new_polytope_venv +*.json +polytope_feature/datacube/quadtree/target +venv_python3_11 +*.txt \ No newline at end of file diff --git a/examples/3D_shipping_route.py b/examples/3D_shipping_route.py index 50db10aa2..571f3a046 100644 --- a/examples/3D_shipping_route.py +++ b/examples/3D_shipping_route.py @@ -6,9 +6,9 @@ import numpy as np from earthkit import data -from polytope.engine.hullslicer import HullSlicer -from polytope.polytope import Polytope, Request -from polytope.shapes import Ellipsoid, Path, Select +from polytope_feature.engine.hullslicer import HullSlicer +from polytope_feature.polytope import Polytope, Request +from polytope_feature.shapes import Ellipsoid, Path, Select class Test: diff --git a/examples/3D_shipping_route_wave_model.py b/examples/3D_shipping_route_wave_model.py index e4edabe9f..81463ba51 100644 --- a/examples/3D_shipping_route_wave_model.py +++ b/examples/3D_shipping_route_wave_model.py @@ -4,10 +4,10 @@ import pytest from eccodes import codes_grib_find_nearest, codes_grib_new_from_file -from polytope.datacube.backends.fdb import FDBDatacube -from polytope.engine.hullslicer import HullSlicer -from polytope.polytope import Polytope, Request -from polytope.shapes import Ellipsoid, Path, Select +from polytope_feature.datacube.backends.fdb import FDBDatacube +from polytope_feature.engine.hullslicer import HullSlicer +from polytope_feature.polytope import Polytope, Request +from polytope_feature.shapes import Ellipsoid, Path, Select from tests.helper_functions import download_test_data diff --git a/examples/4D_flight_path.py b/examples/4D_flight_path.py index d2d1f1eaa..0dc0f8552 100644 --- a/examples/4D_flight_path.py +++ b/examples/4D_flight_path.py @@ -4,10 +4,10 @@ from earthkit import data from PIL import Image -from polytope.datacube.backends.xarray import XArrayDatacube -from polytope.engine.hullslicer import HullSlicer -from polytope.polytope import Polytope, Request -from polytope.shapes import Box, Path, Select +from polytope_feature.datacube.backends.xarray import XArrayDatacube +from polytope_feature.engine.hullslicer import HullSlicer +from polytope_feature.polytope import Polytope, Request +from polytope_feature.shapes import Box, Path, Select class Test: diff --git a/examples/country_slicing.py b/examples/country_slicing.py index 4a5ed30db..fa88b3a5f 100644 --- a/examples/country_slicing.py +++ b/examples/country_slicing.py @@ -4,10 +4,10 @@ from earthkit import data from shapely.geometry import shape -from polytope.datacube.backends.xarray import XArrayDatacube -from polytope.engine.hullslicer import HullSlicer -from polytope.polytope import Polytope, Request -from polytope.shapes import Polygon, Select, Union +from polytope_feature.datacube.backends.xarray import XArrayDatacube +from polytope_feature.engine.hullslicer import HullSlicer +from polytope_feature.polytope import Polytope, Request +from polytope_feature.shapes import Polygon, Select, Union class Test: diff --git a/examples/cyclic_route_around_earth.py b/examples/cyclic_route_around_earth.py index 4e971d677..26cf5ef81 100644 --- a/examples/cyclic_route_around_earth.py +++ b/examples/cyclic_route_around_earth.py @@ -3,10 +3,10 @@ import numpy as np from earthkit import data -from polytope.datacube.backends.xarray import XArrayDatacube -from polytope.engine.hullslicer import HullSlicer -from polytope.polytope import Polytope, Request -from polytope.shapes import Box, PathSegment, Select +from polytope_feature.datacube.backends.xarray import XArrayDatacube +from polytope_feature.engine.hullslicer import HullSlicer +from polytope_feature.polytope import Polytope, Request +from polytope_feature.shapes import Box, PathSegment, Select class Test: diff --git a/examples/healpix_grid_box_example.py b/examples/healpix_grid_box_example.py index 6cfbb5e78..4e759e28f 100644 --- a/examples/healpix_grid_box_example.py +++ b/examples/healpix_grid_box_example.py @@ -3,10 +3,10 @@ from earthkit import data from eccodes import codes_grib_find_nearest, codes_grib_new_from_file -from polytope.datacube.backends.xarray import XArrayDatacube -from polytope.engine.hullslicer import HullSlicer -from polytope.polytope import Polytope, Request -from polytope.shapes import Box, Select +from polytope_feature.datacube.backends.xarray import XArrayDatacube +from polytope_feature.engine.hullslicer import HullSlicer +from polytope_feature.polytope import Polytope, Request +from polytope_feature.shapes import Box, Select class TestOctahedralGrid: diff --git a/examples/octahedral_grid_box_example.py b/examples/octahedral_grid_box_example.py index a4cddcacd..8293c9eb7 100644 --- a/examples/octahedral_grid_box_example.py +++ b/examples/octahedral_grid_box_example.py @@ -6,10 +6,10 @@ from eccodes import codes_grib_find_nearest, codes_grib_new_from_file from matplotlib import markers -from polytope.datacube.backends.xarray import XArrayDatacube -from polytope.engine.hullslicer import HullSlicer -from polytope.polytope import Polytope, Request -from polytope.shapes import Box, Select +from polytope_feature.datacube.backends.xarray import XArrayDatacube +from polytope_feature.engine.hullslicer import HullSlicer +from polytope_feature.polytope import Polytope, Request +from polytope_feature.shapes import Box, Select def find_nearest_latlon(grib_file, target_lat, target_lon): diff --git a/examples/octahedral_grid_country_example.py b/examples/octahedral_grid_country_example.py index 280701d2b..9d777872b 100644 --- a/examples/octahedral_grid_country_example.py +++ b/examples/octahedral_grid_country_example.py @@ -7,10 +7,10 @@ from matplotlib import markers from shapely.geometry import shape -from polytope.datacube.backends.xarray import XArrayDatacube -from polytope.engine.hullslicer import HullSlicer -from polytope.polytope import Polytope, Request -from polytope.shapes import Polygon, Select, Union +from polytope_feature.datacube.backends.xarray import XArrayDatacube +from polytope_feature.engine.hullslicer import HullSlicer +from polytope_feature.polytope import Polytope, Request +from polytope_feature.shapes import Polygon, Select, Union def find_nearest_latlon(grib_file, target_lat, target_lon): diff --git a/examples/read_me_example.py b/examples/read_me_example.py index 0c13e1274..41e3247e6 100644 --- a/examples/read_me_example.py +++ b/examples/read_me_example.py @@ -1,8 +1,8 @@ import numpy as np from earthkit import data -from polytope.polytope import Polytope, Request -from polytope.shapes import Box, Select +from polytope_feature.polytope import Polytope, Request +from polytope_feature.shapes import Box, Select ds = data.from_source("file", "./examples/data/winds.grib") array = ds.to_xarray() diff --git a/examples/slicing_all_ecmwf_countries.py b/examples/slicing_all_ecmwf_countries.py index 7178305a9..dd5e029af 100644 --- a/examples/slicing_all_ecmwf_countries.py +++ b/examples/slicing_all_ecmwf_countries.py @@ -4,10 +4,10 @@ from earthkit import data from shapely.geometry import shape -from polytope.datacube.backends.xarray import XArrayDatacube -from polytope.engine.hullslicer import HullSlicer -from polytope.polytope import Polytope, Request -from polytope.shapes import Polygon, Select, Union +from polytope_feature.datacube.backends.xarray import XArrayDatacube +from polytope_feature.engine.hullslicer import HullSlicer +from polytope_feature.polytope import Polytope, Request +from polytope_feature.shapes import Polygon, Select, Union class Test: diff --git a/examples/timeseries_example.py b/examples/timeseries_example.py index bd350cc13..edbed2efc 100644 --- a/examples/timeseries_example.py +++ b/examples/timeseries_example.py @@ -3,10 +3,10 @@ from earthkit import data from shapely.geometry import shape -from polytope.datacube.backends.xarray import XArrayDatacube -from polytope.engine.hullslicer import HullSlicer -from polytope.polytope import Polytope, Request -from polytope.shapes import Polygon, Select, Union +from polytope_feature.datacube.backends.xarray import XArrayDatacube +from polytope_feature.engine.hullslicer import HullSlicer +from polytope_feature.polytope import Polytope, Request +from polytope_feature.shapes import Polygon, Select, Union class Test: diff --git a/examples/wind_farms.py b/examples/wind_farms.py index 23699eb36..7c8284c8a 100644 --- a/examples/wind_farms.py +++ b/examples/wind_farms.py @@ -6,10 +6,10 @@ from osgeo import gdal from shapely.geometry import shape -from polytope.datacube.backends.xarray import XArrayDatacube -from polytope.engine.hullslicer import HullSlicer -from polytope.polytope import Polytope, Request -from polytope.shapes import Polygon, Select, Union +from polytope_feature.datacube.backends.xarray import XArrayDatacube +from polytope_feature.engine.hullslicer import HullSlicer +from polytope_feature.polytope import Polytope, Request +from polytope_feature.shapes import Polygon, Select, Union class Test: diff --git a/performance/fdb_performance.py b/performance/fdb_performance.py index 78819d462..a8d9e0133 100644 --- a/performance/fdb_performance.py +++ b/performance/fdb_performance.py @@ -2,10 +2,10 @@ import pandas as pd -from polytope.datacube.backends.fdb import FDBDatacube -from polytope.engine.hullslicer import HullSlicer -from polytope.polytope import Polytope, Request -from polytope.shapes import Box, Select +from polytope_feature.datacube.backends.fdb import FDBDatacube +from polytope_feature.engine.hullslicer import HullSlicer +from polytope_feature.polytope import Polytope, Request +from polytope_feature.shapes import Box, Select class TestSlicingFDBDatacube: diff --git a/performance/fdb_performance_3D.py b/performance/fdb_performance_3D.py index 88fb85c92..75c87ba62 100644 --- a/performance/fdb_performance_3D.py +++ b/performance/fdb_performance_3D.py @@ -2,10 +2,10 @@ import pandas as pd -from polytope.datacube.backends.fdb import FDBDatacube -from polytope.engine.hullslicer import HullSlicer -from polytope.polytope import Polytope, Request -from polytope.shapes import Box, Select, Span +from polytope_feature.datacube.backends.fdb import FDBDatacube +from polytope_feature.engine.hullslicer import HullSlicer +from polytope_feature.polytope import Polytope, Request +from polytope_feature.shapes import Box, Select, Span class TestSlicingFDBDatacube: diff --git a/performance/fdb_slice_many_numbers_timeseries.py b/performance/fdb_slice_many_numbers_timeseries.py index 6d62df54c..662b5cc38 100644 --- a/performance/fdb_slice_many_numbers_timeseries.py +++ b/performance/fdb_slice_many_numbers_timeseries.py @@ -3,13 +3,12 @@ import pandas as pd import pygribjump as gj -from polytope.polytope import Polytope, Request -from polytope.shapes import All, Point, Select +from polytope_feature.polytope import Polytope, Request +from polytope_feature.shapes import All, Point, Select time1 = time.time() # Create a dataarray with 3 labelled axes using different index types -# config = {"class": "od", "expver": "0001", "levtype": "sfc", "type": "pf"} options = { "axis_config": [ {"axis_name": "step", "transformations": [{"name": "type_change", "type": "int"}]}, @@ -62,7 +61,6 @@ Select("class", ["od"]), Select("stream", ["enfo"]), Select("type", ["pf"]), - # Select("latitude", [0.035149384216], method="surrounding"), Point(["latitude", "longitude"], [[0.04, 0]], method="surrounding"), All("number"), ) diff --git a/performance/performance_many_num_steps.py b/performance/performance_many_num_steps.py index eb47b0a06..0cdcefcb2 100644 --- a/performance/performance_many_num_steps.py +++ b/performance/performance_many_num_steps.py @@ -4,9 +4,9 @@ import pandas as pd from earthkit import data -from polytope.datacube.backends.fdb import FDBDatacube -from polytope.polytope import Polytope, Request -from polytope.shapes import All, Point, Select, Span +from polytope_feature.datacube.backends.fdb import FDBDatacube +from polytope_feature.polytope import Polytope, Request +from polytope_feature.shapes import All, Point, Select, Span time1 = time.time() # Create a dataarray with 3 labelled axes using different index types diff --git a/performance/scalability_test.py b/performance/scalability_test.py index ab3c5b299..c0aba173f 100644 --- a/performance/scalability_test.py +++ b/performance/scalability_test.py @@ -3,10 +3,10 @@ import numpy as np import xarray as xr -from polytope.datacube.backends.xarray import XArrayDatacube -from polytope.engine.hullslicer import HullSlicer -from polytope.polytope import Polytope, Request -from polytope.shapes import Box, Disk, Ellipsoid, Select +from polytope_feature.datacube.backends.xarray import XArrayDatacube +from polytope_feature.engine.hullslicer import HullSlicer +from polytope_feature.polytope import Polytope, Request +from polytope_feature.shapes import Box, Disk, Ellipsoid, Select class Test: diff --git a/performance/scalability_test_2.py b/performance/scalability_test_2.py index b10e86f4c..7ba0bff4a 100644 --- a/performance/scalability_test_2.py +++ b/performance/scalability_test_2.py @@ -3,10 +3,10 @@ import numpy as np import xarray as xr -from polytope.datacube.backends.xarray import XArrayDatacube -from polytope.engine.hullslicer import HullSlicer -from polytope.polytope import Polytope, Request -from polytope.shapes import Box, Select, Union +from polytope_feature.datacube.backends.xarray import XArrayDatacube +from polytope_feature.engine.hullslicer import HullSlicer +from polytope_feature.polytope import Polytope, Request +from polytope_feature.shapes import Box, Select, Union class Test: diff --git a/performance_unstructured/octahedral_vs_unstructured_slicing.py b/performance_unstructured/octahedral_vs_unstructured_slicing.py new file mode 100644 index 000000000..ab339f95b --- /dev/null +++ b/performance_unstructured/octahedral_vs_unstructured_slicing.py @@ -0,0 +1,356 @@ +# Compare octahedral grid treated as an unstructured vs as an iso-latitude grid + + +import math +import time + +import pandas as pd +import pygribjump as gj +import xarray as xr + +from polytope_feature.polytope import Polytope, Request +from polytope_feature.shapes import Box, Select, Polygon, ConvexPolytope + +import os + +from earthkit import data + + +import os + +import requests +from eccodes import codes_grib_find_nearest, codes_grib_new_from_file + + +# import logging +# logger = logging.getLogger('') +# logger.setLevel(logging.DEBUG) + + +class HTTPError(Exception): + def __init__(self, status_code, message): + self.status_code = status_code + self.message = message + super().__init__(f"HTTPError {status_code}: {message}") + + +def download_test_data(nexus_url, filename): + local_directory = "./tests/data" + + if not os.path.exists(local_directory): + os.makedirs(local_directory) + + # Construct the full path for the local file + local_file_path = os.path.join(local_directory, filename) + + if not os.path.exists(local_file_path): + session = requests.Session() + response = session.get(nexus_url) + if response.status_code != 200: + raise HTTPError(response.status_code, "Failed to download data.") + # Save the downloaded data to the local file + with open(local_file_path, "wb") as f: + f.write(response.content) + + +# from tests.helper_functions import download_test_data +# os.environ["FDB_HOME"] = "/Users/male/git/fdb-new-home" +# slicer_type = "quadtree" +# file_name = "../../Downloads/icon_grid_0026_R03B07_G.nc" +nexus_url = "https://get.ecmwf.int/test-data/polytope/test-data/foo.grib" +download_test_data(nexus_url, "foo.grib") + +file_name = "/Users/male/git/polytope/tests/data/foo.grib" + +ds = data.from_source("file", "./tests/data/foo.grib") + + +def get_engine_options(slicer_type): + engine_options = { + "step": "hullslicer", + "date": "hullslicer", + "levtype": "hullslicer", + "param": "hullslicer", + "class": "hullslicer", + "domain": "hullslicer", + "expver": "hullslicer", + "stream": "hullslicer", + "type": "hullslicer", + "latitude": slicer_type, + "longitude": slicer_type, + } + + return engine_options + + +def get_grid_points(file_name): + # grid = xr.open_dataset(file_name, engine="netcdf4") + grid = data.from_source("file", file_name).to_xarray() + + longitudes = grid.longitude.values + + # print(longitudes) + latitudes = grid.latitude.values + + points = list(zip(latitudes, longitudes)) + # print(points) + return points + + +def set_up_slicing_unstructured(file_name): + options = { + "axis_config": [ + {"axis_name": "step", "transformations": [{"name": "type_change", "type": "int"}]}, + { + "axis_name": "date", + "transformations": [{"name": "merge", "other_axis": "time", "linkers": ["T", "00"]}], + }, + { + "axis_name": "values", + "transformations": [ + {"name": "mapper", "type": "irregular", "resolution": 0, "axes": [ + "latitude", "longitude"], "md5_hash": "158db321ae8e773681eeb40e0a3d350f"} + ], + }, + {"axis_name": "latitude", "transformations": [{"name": "reverse", "is_reverse": True}]}, + {"axis_name": "longitude", "transformations": [{"name": "cyclic", "range": [0, 360]}]}, + ], + "compressed_axes_config": [ + "longitude", + "latitude", + "levtype", + "date", + "domain", + "expver", + "param", + "class", + "stream", + "type", + ], + "pre_path": {"class": "od", "expver": "0001", "levtype": "sfc", "type": "fc", "stream": "oper"}, + } + + fdbdatacube = gj.GribJump() + + engine_options = get_engine_options("quadtree") + points = get_grid_points(file_name) + + API = Polytope( + datacube=fdbdatacube, + options=options, + engine_options=engine_options, + point_cloud_options=points, + ) + + return API + + +def set_up_slicing_structured(): + options = { + # "axis_config": [ + # {"axis_name": "step", "transformations": [{"name": "type_change", "type": "int"}]}, + # { + # "axis_name": "date", + # "transformations": [{"name": "merge", "other_axis": "time", "linkers": ["T", "00"]}], + # }, + # { + # "axis_name": "values", + # "transformations": [ + # {"name": "mapper", "type": "octahedral", "resolution": 1280, "axes": [ + # "latitude", "longitude"]} + # ], + # }, + # ], + # "pre_path": {"date": "20250110"}, + "axis_config": [ + {"axis_name": "step", "transformations": [{"name": "type_change", "type": "int"}]}, + { + "axis_name": "date", + "transformations": [{"name": "merge", "other_axis": "time", "linkers": ["T", "00"]}], + }, + { + "axis_name": "values", + "transformations": [ + {"name": "mapper", "type": "octahedral", "resolution": 1280, "axes": ["latitude", "longitude"]} + ], + }, + {"axis_name": "latitude", "transformations": [{"name": "reverse", "is_reverse": True}]}, + {"axis_name": "longitude", "transformations": [{"name": "cyclic", "range": [0, 360]}]}, + ], + "compressed_axes_config": [ + "longitude", + "latitude", + "levtype", + "date", + "domain", + "expver", + "param", + "class", + "stream", + "type", + ], + "pre_path": {"class": "od", "expver": "0001", "levtype": "sfc", "type": "fc", "stream": "oper"}, + } + + fdbdatacube = gj.GribJump() + + engine_options = get_engine_options("hullslicer") + + API = Polytope( + datacube=fdbdatacube, + options=options, + engine_options=engine_options, + ) + + return API + + +def request(box_size): + + request = Request( + Select("step", [0]), + Select("levtype", ["sfc"]), + Select("date", [pd.Timestamp("20240103T0000")]), + Select("domain", ["g"]), + Select("expver", ["0001"]), + Select("param", ["167"]), + Select("class", ["od"]), + Select("stream", ["oper"]), + Select("type", ["fc"]), + Box(["latitude", "longitude"], [0, 0], [box_size, box_size]), + # ConvexPolytope(["latitude", "longitude"], [[0, 0], [0, box_size], [box_size, box_size], [box_size, 0]]) + ) + + return request + + +# Compare boxes of different sizes + +box_size_max = 80 + +API_structured = set_up_slicing_structured() +API_unstructured = set_up_slicing_unstructured(file_name) + +request_ = request(box_size_max) + +time1 = time.time() +print("\n\n") + +result = API_structured.retrieve(request_) +print("##################################################") + +print("TIME TO RETRIEVE WITH STRUCTURED SLICER") +print("##################################################") + +print("TIME TO RETRIEVE") +print(time.time() - time1) + +print("NUM LEAVES") +leaves = result.leaves +print(len(leaves)) +print("NUM FOUND VALUES") +tot_num = 0 +for leaf in leaves: + tot_num += len(leaf.result) +print(tot_num) + +time2 = time.time() + +print("\n\n") +result = API_unstructured.retrieve(request_) +print("##################################################") +print("TIME TO RETRIEVE WITH UNSTRUCTURED SLICER") +print("##################################################") +print("TIME TO RETRIEVE") +print(time.time() - time2) + +print("NUM LEAVES") +print(len(result.leaves)) + + +# # Compare polygons with n vertices + +# def find_circle_points(num_points, centre, radius): + +# expanded_radius = _expansion_to_circumscribe_circle(num_points) +# points = _points_on_circle(num_points, expanded_radius) + +# for i in range(0, len(points)): +# x = centre[0] + points[i][0] * radius[0] +# y = centre[1] + points[i][1] * radius[1] +# points[i] = [x, y] +# return points + + +# def _points_on_circle(n, r): +# return [[math.cos(2 * math.pi / n * x) * r, math.sin(2 * math.pi / n * x) * r] for x in range(0, n)] + + +# def _expansion_to_circumscribe_circle(n): +# half_angle_between_segments = math.pi / n +# return 1 / math.cos(half_angle_between_segments) + + +# def request_disk(num_sides): + +# points_disk = find_circle_points(num_sides, [12, 12], [10, 10]) + +# request = Request( +# Select("step", [0]), +# Select("levtype", ["sfc"]), +# Select("date", [pd.Timestamp("20240103T0000")]), +# Select("domain", ["g"]), +# Select("expver", ["0001"]), +# Select("param", ["167"]), +# Select("class", ["od"]), +# Select("stream", ["oper"]), +# Select("type", ["fc"]), +# # Box(["latitude", "longitude"], [0, 0], box_size), +# ConvexPolytope(["latitude", "longitude"], points_disk) +# ) + +# return request + + +# n_sides = 64 + +# API_structured = set_up_slicing_structured() +# API_unstructured = set_up_slicing_unstructured(file_name) + +# request_ = request_disk(n_sides) + +# time1 = time.time() +# print("\n\n") + +# result = API_structured.retrieve(request_) +# print("##################################################") + +# print("TIME TO RETRIEVE WITH STRUCTURED SLICER") +# print("##################################################") + +# print("TIME TO RETRIEVE") +# print(time.time() - time1) + +# print("NUM LEAVES") +# leaves = result.leaves +# print(len(leaves)) +# print("NUM FOUND VALUES") +# tot_num = 0 +# for leaf in leaves: +# tot_num += len(leaf.result) +# print(tot_num) +# result.pprint() + +# time2 = time.time() + +# print("\n\n") +# result = API_unstructured.retrieve(request_) +# print("##################################################") +# print("TIME TO RETRIEVE WITH UNSTRUCTURED SLICER") +# print("##################################################") +# print("TIME TO RETRIEVE") +# print(time.time() - time2) + +# print("NUM LEAVES") +# print(len(result.leaves)) +# result.pprint() diff --git a/performance_unstructured/plot_structured_vs_unstructured_slicing.py b/performance_unstructured/plot_structured_vs_unstructured_slicing.py new file mode 100644 index 000000000..a2af0abd8 --- /dev/null +++ b/performance_unstructured/plot_structured_vs_unstructured_slicing.py @@ -0,0 +1,33 @@ +import matplotlib.pyplot as plt + + +box_size = [10, 20, 40, 80] + +quadtree_build_time = [10.954054832458496, 10.95795202255249, 10.915708065032959, 11.012521028518677] + +num_vals_found = [19226, 72485, 253205, 724400] + +gj_retrieve_structured_time = [0.0176, 0.021518, 0.022053, 0.051661] + +gj_retrieve_unstructured_time = [0.097058, 1, 13, 128] + +tot_retrieve_time_structured = [0.2735710144042969, 0.5628271102905273, 1.5192251205444336, 4.548176288604736] + +tot_retrieve_time_unstructured = [0.9918620586395264, 5.218957185745239, 33.4666748046875, 210.05072903633118] + + +pure_poly_time_structured = [tot_retrieve_time_structured[i] - gj_retrieve_structured_time[i] for i in range(4)] +pure_poly_time_unstructured = [tot_retrieve_time_unstructured[i] - gj_retrieve_unstructured_time[i] for i in range(4)] + +fig, ax = plt.subplots(1, 1) + +# ax.plot(box_size, tot_retrieve_time_structured, color="b") +# ax.plot(box_size, tot_retrieve_time_unstructured, color="r") +# ax.plot(box_size, gj_retrieve_structured_time, color="b", linestyle="--") +# ax.plot(box_size, gj_retrieve_unstructured_time, color="r", linestyle="--") +# ax.set_yscale("log") + +ax.plot(box_size, pure_poly_time_structured, color="b") +ax.plot(box_size, pure_poly_time_unstructured, color="r") + +plt.show() diff --git a/polytope_feature/datacube/backends/datacube.py b/polytope_feature/datacube/backends/datacube.py index 97a8c7c77..9f961b174 100644 --- a/polytope_feature/datacube/backends/datacube.py +++ b/polytope_feature/datacube/backends/datacube.py @@ -152,20 +152,40 @@ def remap_path(self, path: DatacubePath): return path @staticmethod - def create(datacube, config={}, axis_options={}, compressed_axes_options=[], alternative_axes=[], context=None): + def create( + datacube, + config={}, + axis_options={}, + compressed_axes_options=[], + point_cloud_options=None, + alternative_axes=[], + context=None, + ): # TODO: get the configs as None for pre-determined value and change them to empty dictionary inside the function if type(datacube).__name__ == "DataArray": from .xarray import XArrayDatacube - xadatacube = XArrayDatacube(datacube, axis_options, compressed_axes_options, context) + xadatacube = XArrayDatacube(datacube, axis_options, compressed_axes_options, point_cloud_options, context) return xadatacube if type(datacube).__name__ == "GribJump": from .fdb import FDBDatacube fdbdatacube = FDBDatacube( - datacube, config, axis_options, compressed_axes_options, alternative_axes, context + datacube, + config, + axis_options, + compressed_axes_options, + point_cloud_options, + alternative_axes, + context, ) return fdbdatacube + if type(datacube).__name__ == "MockDatacube": + return datacube def check_branching_axes(self, request): pass + + @abstractmethod + def find_point_cloud(self): + pass diff --git a/polytope_feature/datacube/backends/fdb.py b/polytope_feature/datacube/backends/fdb.py index bd60665af..4b53d0634 100644 --- a/polytope_feature/datacube/backends/fdb.py +++ b/polytope_feature/datacube/backends/fdb.py @@ -10,7 +10,14 @@ class FDBDatacube(Datacube): def __init__( - self, gj, config=None, axis_options=None, compressed_axes_options=[], alternative_axes=[], context=None + self, + gj, + config=None, + axis_options=None, + compressed_axes_options=[], + point_cloud_options=None, + alternative_axes=[], + context=None, ): if config is None: config = {} @@ -23,6 +30,7 @@ def __init__( self.unwanted_path = {} self.axis_options = axis_options + self.has_point_cloud = point_cloud_options # NOTE: here, will be True/False partial_request = config # Find values in the level 3 FDB datacube @@ -70,7 +78,12 @@ def __init__( val = self._axes[name].type self._check_and_add_axes(options, name, val) - logging.info("Polytope created axes for %s", self._axes.keys()) + logging.info("Polytope created axes for: " + str(self._axes.keys())) + + def find_point_cloud(self): + # TODO: somehow, find the point cloud of irregular grid if it exists + if self.has_point_cloud: + return self.has_point_cloud def check_branching_axes(self, request): polytopes = request.polytopes() @@ -167,7 +180,25 @@ def get_fdb_requests( key_value_path, leaf_path, self.unwanted_path ) leaf_path.update(key_value_path) - if len(requests.children[0].children[0].children) == 0: + if not requests.children[0].children: + # TODO: the children are now GridNodes + print("ARE HERE WHERE WE EXPECT??") + print(requests.children[0].indexes) + (path, current_start_idxs, fdb_node_ranges, lat_length) = self.get_leaf_grid_node_values(requests, leaf_path) + print("AND GOT HERE NOW??") + ( + original_indices, + sorted_request_ranges, + fdb_node_ranges, + ) = self.sort_fdb_request_ranges(current_start_idxs, lat_length, fdb_node_ranges, True) + fdb_requests.append((path, sorted_request_ranges)) + fdb_requests_decoding_info.append((original_indices, fdb_node_ranges)) + pass + elif not requests.children[0].children[0].children: + # TODO: We have a GridNode instead of a normal node as a grandchild + for c in requests.children: + self.get_fdb_requests(c, fdb_requests, fdb_requests_decoding_info, leaf_path) + elif len(requests.children[0].children[0].children) == 0: # find the fdb_requests and associated nodes to which to add results (path, current_start_idxs, fdb_node_ranges, lat_length) = self.get_2nd_last_values(requests, leaf_path) ( @@ -253,6 +284,31 @@ def nearest_lat_lon_search(self, requests): if value not in possible_lons: lon_child.remove_compressed_branch(value) + def get_leaf_grid_node_values(self, requests, leaf_path=None): + # here we only work with GridNodes so actually it's a completely "flat" leaf + if leaf_path is None: + leaf_path = {} + + # TODO: find the nearest lat/lon + + tree_length = len(requests.children) + current_start_idxs = [False] * tree_length + fdb_node_ranges = [False] * tree_length + + for i, c in enumerate(requests.children): + # Add each child to the fdb_node_ranges + fdb_node_ranges[i] = [[c]] + # TODO: unmap each lat/lon successively with the right datacube axis + current_start_idxs[i] = [c.indexes] + # TODO: update the leaf path + # key_value_path = {c.axes[0]: (c.values[0],), c.axes[1]: (c.values[1],)} + # leaf_path.update(key_value_path) + + # TODO: remove the leaf_path_copy leaves we do not need + # TODO: return everything like before + leaf_path_copy = deepcopy(leaf_path) + return (leaf_path_copy, current_start_idxs, fdb_node_ranges, tree_length) + def get_2nd_last_values(self, requests, leaf_path=None): if leaf_path is None: leaf_path = {} @@ -282,6 +338,7 @@ def get_2nd_last_values(self, requests, leaf_path=None): leaf_path_copy = deepcopy(leaf_path) leaf_path_copy.pop("values", None) + leaf_path_copy.pop("index") return (leaf_path_copy, current_start_idxs, fdb_node_ranges, lat_length) def get_last_layer_before_leaf(self, requests, leaf_path, current_idx, fdb_range_n): @@ -290,6 +347,7 @@ def get_last_layer_before_leaf(self, requests, leaf_path, current_idx, fdb_range for i, c in enumerate(requests.children): # now c are the leaves of the initial tree key_value_path = {c.axis.name: c.values} + leaf_path["index"] = c.indexes ax = c.axis (key_value_path, leaf_path, self.unwanted_path) = ax.unmap_path_key( key_value_path, leaf_path, self.unwanted_path @@ -317,43 +375,80 @@ def assign_fdb_output_to_nodes(self, output_values, fdb_requests_decoding_info): interm_request_output_values = request_output_values[0][i][0] n.result.extend(interm_request_output_values) - def sort_fdb_request_ranges(self, current_start_idx, lat_length, fdb_node_ranges): - (new_fdb_node_ranges, new_current_start_idx) = self.remove_duplicates_in_request_ranges( - fdb_node_ranges, current_start_idx - ) - interm_request_ranges = [] - # TODO: modify the start indexes to have as many arrays as the request ranges - new_fdb_node_ranges = [] - for i in range(lat_length): - interm_fdb_nodes = fdb_node_ranges[i] - old_interm_start_idx = current_start_idx[i] - for j in range(len(old_interm_start_idx)): - # TODO: if we sorted the cyclic values in increasing order on the tree too, - # then we wouldn't have to sort here? - sorted_list = sorted(enumerate(old_interm_start_idx[j]), key=lambda x: x[1]) - original_indices_idx, interm_start_idx = zip(*sorted_list) - for interm_fdb_nodes_obj in interm_fdb_nodes[j]: - interm_fdb_nodes_obj.values = tuple([interm_fdb_nodes_obj.values[k] for k in original_indices_idx]) - if abs(interm_start_idx[-1] + 1 - interm_start_idx[0]) <= len(interm_start_idx): - current_request_ranges = (interm_start_idx[0], interm_start_idx[-1] + 1) - interm_request_ranges.append(current_request_ranges) - new_fdb_node_ranges.append(interm_fdb_nodes[j]) - else: - jumps = list(map(operator.sub, interm_start_idx[1:], interm_start_idx[:-1])) - last_idx = 0 - for k, jump in enumerate(jumps): - if jump > 1: - current_request_ranges = (interm_start_idx[last_idx], interm_start_idx[k] + 1) - new_fdb_node_ranges.append(interm_fdb_nodes[j]) - last_idx = k + 1 - interm_request_ranges.append(current_request_ranges) - if k == len(interm_start_idx) - 2: - current_request_ranges = (interm_start_idx[last_idx], interm_start_idx[-1] + 1) - interm_request_ranges.append(current_request_ranges) - new_fdb_node_ranges.append(interm_fdb_nodes[j]) - request_ranges_with_idx = list(enumerate(interm_request_ranges)) - sorted_list = sorted(request_ranges_with_idx, key=lambda x: x[1][0]) - original_indices, sorted_request_ranges = zip(*sorted_list) + def sort_fdb_request_ranges(self, current_start_idx, lat_length, fdb_node_ranges, is_grid_node=False): + if not is_grid_node: + (new_fdb_node_ranges, new_current_start_idx) = self.remove_duplicates_in_request_ranges( + fdb_node_ranges, current_start_idx + ) + interm_request_ranges = [] + # TODO: modify the start indexes to have as many arrays as the request ranges + new_fdb_node_ranges = [] + for i in range(lat_length): + interm_fdb_nodes = fdb_node_ranges[i] + old_interm_start_idx = current_start_idx[i] + for j in range(len(old_interm_start_idx)): + # TODO: if we sorted the cyclic values in increasing order on the tree too, + # then we wouldn't have to sort here? + sorted_list = sorted(enumerate(old_interm_start_idx[j]), key=lambda x: x[1]) + original_indices_idx, interm_start_idx = zip(*sorted_list) + for interm_fdb_nodes_obj in interm_fdb_nodes[j]: + interm_fdb_nodes_obj.values = tuple([interm_fdb_nodes_obj.values[k] + for k in original_indices_idx]) + if abs(interm_start_idx[-1] + 1 - interm_start_idx[0]) <= len(interm_start_idx): + current_request_ranges = (interm_start_idx[0], interm_start_idx[-1] + 1) + interm_request_ranges.append(current_request_ranges) + new_fdb_node_ranges.append(interm_fdb_nodes[j]) + else: + jumps = list(map(operator.sub, interm_start_idx[1:], interm_start_idx[:-1])) + last_idx = 0 + for k, jump in enumerate(jumps): + if jump > 1: + current_request_ranges = (interm_start_idx[last_idx], interm_start_idx[k] + 1) + new_fdb_node_ranges.append(interm_fdb_nodes[j]) + last_idx = k + 1 + interm_request_ranges.append(current_request_ranges) + if k == len(interm_start_idx) - 2: + current_request_ranges = (interm_start_idx[last_idx], interm_start_idx[-1] + 1) + interm_request_ranges.append(current_request_ranges) + new_fdb_node_ranges.append(interm_fdb_nodes[j]) + request_ranges_with_idx = list(enumerate(interm_request_ranges)) + sorted_list = sorted(request_ranges_with_idx, key=lambda x: x[1][0]) + original_indices, sorted_request_ranges = zip(*sorted_list) + else: + interm_request_ranges = [] + # TODO: modify the start indexes to have as many arrays as the request ranges + new_fdb_node_ranges = [] + for i in range(lat_length): + interm_fdb_nodes = fdb_node_ranges[i] + old_interm_start_idx = current_start_idx[i] + for j in range(len(old_interm_start_idx)): + # TODO: if we sorted the cyclic values in increasing order on the tree too, + # then we wouldn't have to sort here? + sorted_list = sorted(enumerate(old_interm_start_idx[j]), key=lambda x: x[1]) + original_indices_idx, interm_start_idx = zip(*sorted_list) + # for interm_fdb_nodes_obj in interm_fdb_nodes[j]: + # interm_fdb_nodes_obj.values = tuple([interm_fdb_nodes_obj.values[k] + # for k in original_indices_idx]) + if abs(interm_start_idx[-1] + 1 - interm_start_idx[0]) <= len(interm_start_idx): + current_request_ranges = (interm_start_idx[0], interm_start_idx[-1] + 1) + interm_request_ranges.append(current_request_ranges) + new_fdb_node_ranges.append(interm_fdb_nodes[j]) + else: + jumps = list(map(operator.sub, interm_start_idx[1:], interm_start_idx[:-1])) + last_idx = 0 + for k, jump in enumerate(jumps): + if jump > 1: + current_request_ranges = (interm_start_idx[last_idx], interm_start_idx[k] + 1) + new_fdb_node_ranges.append(interm_fdb_nodes[j]) + last_idx = k + 1 + interm_request_ranges.append(current_request_ranges) + if k == len(interm_start_idx) - 2: + current_request_ranges = (interm_start_idx[last_idx], interm_start_idx[-1] + 1) + interm_request_ranges.append(current_request_ranges) + new_fdb_node_ranges.append(interm_fdb_nodes[j]) + request_ranges_with_idx = list(enumerate(interm_request_ranges)) + sorted_list = sorted(request_ranges_with_idx, key=lambda x: x[1][0]) + original_indices, sorted_request_ranges = zip(*sorted_list) return (original_indices, sorted_request_ranges, new_fdb_node_ranges) def datacube_natural_indexes(self, axis, subarray): diff --git a/polytope_feature/datacube/backends/mock.py b/polytope_feature/datacube/backends/mock.py index c23b80802..f652bfe99 100644 --- a/polytope_feature/datacube/backends/mock.py +++ b/polytope_feature/datacube/backends/mock.py @@ -13,10 +13,10 @@ def __init__(self, dimensions, compressed_axes_options=[]): self.dimensions = dimensions - self.mappers = {} + self._axes = {} for name in self.dimensions: - self.mappers[name] = deepcopy(IntDatacubeAxis()) - self.mappers[name].name = name + self._axes[name] = deepcopy(IntDatacubeAxis()) + self._axes[name].name = name self.stride = {} stride_cumulative = 1 @@ -42,7 +42,7 @@ def get(self, requests: TensorIndexTree, context=None): r.remove_branch() def get_mapper(self, axis): - return self.mappers[axis] + return self._axes[axis] def get_indices(self, path: DatacubePath, axis, lower, upper, method=None): if lower == upper == math.ceil(lower): @@ -59,7 +59,7 @@ def has_index(self, path: DatacubePath, axis, index): @property def axes(self): - return self.mappers + return self._axes def validate(self, axes): return validate_axes(self.axes, axes) @@ -69,3 +69,6 @@ def ax_vals(self, name): def _find_indexes_between(self, axis, indexes, low, up): pass + + def find_point_cloud(self): + pass diff --git a/polytope_feature/datacube/backends/xarray.py b/polytope_feature/datacube/backends/xarray.py index 735c87e25..67f3e1cb5 100644 --- a/polytope_feature/datacube/backends/xarray.py +++ b/polytope_feature/datacube/backends/xarray.py @@ -9,7 +9,14 @@ class XArrayDatacube(Datacube): """Xarray arrays are labelled, axes can be defined as strings or integers (e.g. "time" or 0).""" - def __init__(self, dataarray: xr.DataArray, axis_options=None, compressed_axes_options=[], context=None): + def __init__( + self, + dataarray: xr.DataArray, + axis_options=None, + compressed_axes_options=[], + point_cloud_options=None, + context=None, + ): super().__init__(axis_options, compressed_axes_options) if axis_options is None: axis_options = {} @@ -17,6 +24,7 @@ def __init__(self, dataarray: xr.DataArray, axis_options=None, compressed_axes_o self.axis_counter = 0 self._axes = None self.dataarray = dataarray + self.has_point_cloud = point_cloud_options for name, values in dataarray.coords.variables.items(): options = None @@ -50,9 +58,12 @@ def __init__(self, dataarray: xr.DataArray, axis_options=None, compressed_axes_o val = self._axes[name].type self._check_and_add_axes(options, name, val) + def find_point_cloud(self): + # TODO: somehow, find the point cloud of irregular grid if it exists + if self.has_point_cloud: + return self.has_point_cloud + def get(self, requests, context=None, leaf_path=None, axis_counter=0): - if context is None: - context = {} if leaf_path is None: leaf_path = {} if requests.axis.name == "root": @@ -68,6 +79,8 @@ def get(self, requests, context=None, leaf_path=None, axis_counter=0): if len(requests.children) != 0: # We are not a leaf and we loop over for c in requests.children: + if axis_counter == self.axis_counter - 1: + leaf_path["index"] = c.indexes self.get(c, context, leaf_path, axis_counter + 1) else: if self.axis_counter != axis_counter: @@ -78,7 +91,8 @@ def get(self, requests, context=None, leaf_path=None, axis_counter=0): unmapped_path = {} self.refit_path(leaf_path_copy, unmapped_path, leaf_path) for key in leaf_path_copy: - leaf_path_copy[key] = list(leaf_path_copy[key]) + if isinstance(leaf_path_copy[key], tuple): + leaf_path_copy[key] = list(leaf_path_copy[key]) for key in unmapped_path: if isinstance(unmapped_path[key], tuple): unmapped_path[key] = list(unmapped_path[key]) @@ -108,9 +122,9 @@ def refit_path(self, path_copy, unmapped_path, path): for key in path.keys(): if key not in self.dataarray.dims: path_copy.pop(key) - if key not in self.dataarray.coords.dtypes: + elif key not in self.dataarray.coords.dtypes: unmapped_path.update({key: path[key]}) - path_copy.pop(key) + path_copy.pop(key, None) for key in self.dataarray.coords.dtypes: key_dtype = self.dataarray.coords.dtypes[key] if key_dtype.type is np.str_ and key in path.keys(): diff --git a/polytope_feature/datacube/quad_tree.py b/polytope_feature/datacube/quad_tree.py new file mode 100644 index 000000000..96c1ba44b --- /dev/null +++ b/polytope_feature/datacube/quad_tree.py @@ -0,0 +1,176 @@ +from ..engine.hullslicer import slice +from ..engine.slicing_tools import slice_in_two + +""" + + QuadTree data structure + + comes from: https://github.com/mdrasmus/compbio/blob/master/rasmus/quadtree.py + + is specific to 2D + +""" + + +class QuadNode: + def __init__(self, item, index): + self.item = item + self.index = index + + def is_contained_in(self, polygon): + # implement method to check if the node point is inside the polygon + node_x, node_y = self.item + + sliced_vertical_polygon = slice(polygon, polygon._axes[0], node_x, 0) + if sliced_vertical_polygon: + lower_y, upper_y, idx = sliced_vertical_polygon.extents(polygon._axes[1]) + if lower_y <= node_y <= upper_y: + return True + return False + + +class QuadTree: + # TODO: do we need the max_depth? + MAX = 3 + MAX_DEPTH = 20 + + def __init__(self, x=0, y=0, size=[180, 90], depth=0): + self.nodes = [] + self.children = [] + self.center = (x, y) + self.size = tuple(size) + self.depth = depth + self.node_items = set() + + def quadrant_rectangle_points(self): + return set( + [ + (self.center[0] + (self.size[0]), self.center[1] + (self.size[1])), + (self.center[0] + (self.size[0]), self.center[1] - (self.size[1])), + (self.center[0] - (self.size[0]), self.center[1] + (self.size[1])), + (self.center[0] - (self.size[0]), self.center[1] - (self.size[1])), + ] + ) + + def build_point_tree(self, points): + # TODO: SLOW, scales linearly with number of points + for index, p in enumerate(points): + self.insert(tuple(p), index) + + def pprint(self): + if self.depth == 0: + print("\n") + if len(self.children) == 0: + for n in self.nodes: + print("\t" * (self.depth + 1) + "\u21b3" + str(n.item[0]) + " , " + str(n.item[1])) + for child in self.children: + print( + "\t" * (self.depth + 1) + + "\u21b3" + + str(child.center[0] - child.size[0]) + + " , " + + str(child.center[1] - child.size[1]) + + " , " + + str(child.center[0] + child.size[0]) + + " , " + + str(child.center[1] + child.size[1]) + ) + child.pprint() + + def insert(self, item, index): + if not self.children: + node = QuadNode(item, index) + if item not in self.node_items: + self.nodes.append(node) + self.node_items.add(node.item) + + if len(self.nodes) > self.MAX and self.depth < self.MAX_DEPTH: + self.split() + return node + else: + return self.insert_into_children(item, index) + + def insert_into_children(self, item, index): + x, y = item + cx, cy = self.center + # try to insert into children + if x <= cx: + if y <= cy: + self.children[0].insert(item, index) + if y >= cy: + self.children[1].insert(item, index) + if x >= cx: + if y <= cy: + self.children[2].insert(item, index) + if y >= cy: + self.children[3].insert(item, index) + + def split(self): + half_size = [s / 2 for s in self.size] + x_center, y_center = self.center[0], self.center[1] + hx, hy = half_size + + new_centers = [ + (x_center - hx, y_center - hy), + (x_center - hx, y_center + hy), + (x_center + hx, y_center - hy), + (x_center + hx, y_center + hy), + ] + + self.children = [ + QuadTree(new_center[0], new_center[1], half_size, self.depth + 1) for new_center in new_centers + ] + + nodes = self.nodes + self.nodes = [] + for node in nodes: + self.insert_into_children(node.item, node.index) + + def query_polygon(self, polygon, results=None): + # intersect quad tree with a 2D polygon + if results is None: + results = set() + + # intersect the children with the polygon + # TODO: here, we create None polygons... think about how to handle them + if polygon is None: + pass + else: + polygon_points = set([tuple(point) for point in polygon.points]) + # TODO: are these the right points which we are comparing, ie the points on the polygon + # and the points on the rectangle quadrant? + if polygon_points == self.quadrant_rectangle_points(): + for node in self.find_nodes_in(): + results.add(node) + else: + if len(self.children) > 0: + # first slice vertically + left_polygon, right_polygon = slice_in_two(polygon, self.center[0], 0) + + # then slice horizontally + # ie need to slice the left and right polygons each in two to have the 4 quadrant polygons + + q1_polygon, q2_polygon = slice_in_two(left_polygon, self.center[1], 1) + q3_polygon, q4_polygon = slice_in_two(right_polygon, self.center[1], 1) + + # now query these 4 polygons further down the quadtree + self.children[0].query_polygon(q1_polygon, results) + self.children[1].query_polygon(q2_polygon, results) + self.children[2].query_polygon(q3_polygon, results) + self.children[3].query_polygon(q4_polygon, results) + + results.update(node for node in self.nodes if node.is_contained_in(polygon)) + + return results + + def find_nodes_in(self, results=None): + # TODO: find the nodes that are in this subtree + if results is None: + results = set() + if len(self.children) > 0: + # there are children which we need to iterate through + for child in self.children: + child.find_nodes_in(results) + for node in self.nodes: + results.add(node) + return results diff --git a/polytope_feature/datacube/quadtree/Cargo.lock b/polytope_feature/datacube/quadtree/Cargo.lock new file mode 100644 index 000000000..31465fd8d --- /dev/null +++ b/polytope_feature/datacube/quadtree/Cargo.lock @@ -0,0 +1,620 @@ +# This file is automatically @generated by Cargo. +# It is not intended for manual editing. +version = 4 + +[[package]] +name = "aho-corasick" +version = "1.1.3" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "8e60d3430d3a69478ad0993f19238d2df97c507009a52b3c10addcd7f6bcb916" +dependencies = [ + "memchr", +] + +[[package]] +name = "autocfg" +version = "1.4.0" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "ace50bade8e6234aa140d9a2f552bbee1db4d353f69b8217bc503490fc1a9f26" + +[[package]] +name = "bindgen" +version = "0.69.5" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "271383c67ccabffb7381723dea0672a673f292304fcb45c01cc648c7a8d58088" +dependencies = [ + "bitflags", + "cexpr", + "clang-sys", + "itertools", + "lazy_static", + "lazycell", + "log", + "prettyplease", + "proc-macro2", + "quote", + "regex", + "rustc-hash", + "shlex", + "syn", + "which", +] + +[[package]] +name = "bitflags" +version = "2.8.0" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "8f68f53c83ab957f72c32642f3868eec03eb974d1fb82e453128456482613d36" + +[[package]] +name = "cc" +version = "1.2.16" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "be714c154be609ec7f5dad223a33bf1482fff90472de28f7362806e6d4832b8c" +dependencies = [ + "shlex", +] + +[[package]] +name = "cexpr" +version = "0.6.0" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "6fac387a98bb7c37292057cffc56d62ecb629900026402633ae9160df93a8766" +dependencies = [ + "nom", +] + +[[package]] +name = "cfg-if" +version = "1.0.0" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "baf1de4339761588bc0619e3cbc0120ee582ebb74b53b4efbf79117bd2da40fd" + +[[package]] +name = "clang-sys" +version = "1.8.1" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "0b023947811758c97c59bf9d1c188fd619ad4718dcaa767947df1cadb14f39f4" +dependencies = [ + "glob", + "libc", + "libloading", +] + +[[package]] +name = "crossbeam-deque" +version = "0.8.6" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "9dd111b7b7f7d55b72c0a6ae361660ee5853c9af73f70c3c2ef6858b950e2e51" +dependencies = [ + "crossbeam-epoch", + "crossbeam-utils", +] + +[[package]] +name = "crossbeam-epoch" +version = "0.9.18" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "5b82ac4a3c2ca9c3460964f020e1402edd5753411d7737aa39c3714ad1b5420e" +dependencies = [ + "crossbeam-utils", +] + +[[package]] +name = "crossbeam-utils" +version = "0.8.21" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "d0a5c400df2834b80a4c3327b3aad3a4c4cd4de0629063962b03235697506a28" + +[[package]] +name = "either" +version = "1.15.0" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "48c757948c5ede0e46177b7add2e67155f70e33c07fea8284df6576da70b3719" + +[[package]] +name = "errno" +version = "0.3.10" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "33d852cb9b869c2a9b3df2f71a3074817f01e1844f839a144f5fcef059a4eb5d" +dependencies = [ + "libc", + "windows-sys", +] + +[[package]] +name = "glob" +version = "0.3.2" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "a8d1add55171497b4705a648c6b583acafb01d58050a51727785f0b2c8e0a2b2" + +[[package]] +name = "heck" +version = "0.4.1" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "95505c38b4572b2d910cecb0281560f54b440a19336cbbcb27bf6ce6adc6f5a8" + +[[package]] +name = "home" +version = "0.5.11" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "589533453244b0995c858700322199b2becb13b627df2851f64a2775d024abcf" +dependencies = [ + "windows-sys", +] + +[[package]] +name = "indoc" +version = "2.0.5" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "b248f5224d1d606005e02c97f5aa4e88eeb230488bcc03bc9ca4d7991399f2b5" + +[[package]] +name = "itertools" +version = "0.12.1" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "ba291022dbbd398a455acf126c1e341954079855bc60dfdda641363bd6922569" +dependencies = [ + "either", +] + +[[package]] +name = "lazy_static" +version = "1.5.0" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "bbd2bcb4c963f2ddae06a2efc7e9f3591312473c50c6685e1f298068316e66fe" + +[[package]] +name = "lazycell" +version = "1.3.0" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "830d08ce1d1d941e6b30645f1a0eb5643013d835ce3779a5fc208261dbe10f55" + +[[package]] +name = "libc" +version = "0.2.169" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "b5aba8db14291edd000dfcc4d620c7ebfb122c613afb886ca8803fa4e128a20a" + +[[package]] +name = "libloading" +version = "0.8.6" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "fc2f4eb4bc735547cfed7c0a4922cbd04a4655978c09b54f1f7b228750664c34" +dependencies = [ + "cfg-if", + "windows-targets", +] + +[[package]] +name = "linux-raw-sys" +version = "0.4.15" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "d26c52dbd32dccf2d10cac7725f8eae5296885fb5703b261f7d0a0739ec807ab" + +[[package]] +name = "lock_api" +version = "0.4.12" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "07af8b9cdd281b7915f413fa73f29ebd5d55d0d3f0155584dade1ff18cea1b17" +dependencies = [ + "autocfg", + "scopeguard", +] + +[[package]] +name = "log" +version = "0.4.26" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "30bde2b3dc3671ae49d8e2e9f044c7c005836e7a023ee57cffa25ab82764bb9e" + +[[package]] +name = "memchr" +version = "2.7.4" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "78ca9ab1a0babb1e7d5695e3530886289c18cf2f87ec19a575a0abdce112e3a3" + +[[package]] +name = "memoffset" +version = "0.9.1" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "488016bfae457b036d996092f6cb448677611ce4449e970ceaf42695203f218a" +dependencies = [ + "autocfg", +] + +[[package]] +name = "minimal-lexical" +version = "0.2.1" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "68354c5c6bd36d73ff3feceb05efa59b6acb7626617f4962be322a825e61f79a" + +[[package]] +name = "nom" +version = "7.1.3" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "d273983c5a657a70a3e8f2a01329822f3b8c8172b73826411a55751e404a0a4a" +dependencies = [ + "memchr", + "minimal-lexical", +] + +[[package]] +name = "once_cell" +version = "1.20.3" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "945462a4b81e43c4e3ba96bd7b49d834c6f61198356aa858733bc4acf3cbe62e" + +[[package]] +name = "parking_lot" +version = "0.12.3" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "f1bf18183cf54e8d6059647fc3063646a1801cf30896933ec2311622cc4b9a27" +dependencies = [ + "lock_api", + "parking_lot_core", +] + +[[package]] +name = "parking_lot_core" +version = "0.9.10" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "1e401f977ab385c9e4e3ab30627d6f26d00e2c73eef317493c4ec6d468726cf8" +dependencies = [ + "cfg-if", + "libc", + "redox_syscall", + "smallvec", + "windows-targets", +] + +[[package]] +name = "portable-atomic" +version = "1.10.0" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "280dc24453071f1b63954171985a0b0d30058d287960968b9b2aca264c8d4ee6" + +[[package]] +name = "prettyplease" +version = "0.2.30" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "f1ccf34da56fc294e7d4ccf69a85992b7dfb826b7cf57bac6a70bba3494cc08a" +dependencies = [ + "proc-macro2", + "syn", +] + +[[package]] +name = "proc-macro2" +version = "1.0.93" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "60946a68e5f9d28b0dc1c21bb8a97ee7d018a8b322fa57838ba31cc878e22d99" +dependencies = [ + "unicode-ident", +] + +[[package]] +name = "pyo3" +version = "0.20.3" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "53bdbb96d49157e65d45cc287af5f32ffadd5f4761438b527b055fb0d4bb8233" +dependencies = [ + "cfg-if", + "indoc", + "libc", + "memoffset", + "parking_lot", + "portable-atomic", + "pyo3-build-config", + "pyo3-ffi", + "pyo3-macros", + "unindent", +] + +[[package]] +name = "pyo3-build-config" +version = "0.20.3" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "deaa5745de3f5231ce10517a1f5dd97d53e5a2fd77aa6b5842292085831d48d7" +dependencies = [ + "once_cell", + "target-lexicon", +] + +[[package]] +name = "pyo3-ffi" +version = "0.20.3" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "62b42531d03e08d4ef1f6e85a2ed422eb678b8cd62b762e53891c05faf0d4afa" +dependencies = [ + "libc", + "pyo3-build-config", +] + +[[package]] +name = "pyo3-macros" +version = "0.20.3" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "7305c720fa01b8055ec95e484a6eca7a83c841267f0dd5280f0c8b8551d2c158" +dependencies = [ + "proc-macro2", + "pyo3-macros-backend", + "quote", + "syn", +] + +[[package]] +name = "pyo3-macros-backend" +version = "0.20.3" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "7c7e9b68bb9c3149c5b0cade5d07f953d6d125eb4337723c4ccdb665f1f96185" +dependencies = [ + "heck", + "proc-macro2", + "pyo3-build-config", + "quote", + "syn", +] + +[[package]] +name = "qhull" +version = "0.4.0" +dependencies = [ + "qhull-sys", + "thiserror", +] + +[[package]] +name = "qhull-sys" +version = "0.4.0" +dependencies = [ + "bindgen", + "cc", +] + +[[package]] +name = "quadtree" +version = "0.1.0" +dependencies = [ + "pyo3", + "qhull", + "rayon", +] + +[[package]] +name = "quote" +version = "1.0.38" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "0e4dccaaaf89514f546c693ddc140f729f958c247918a13380cccc6078391acc" +dependencies = [ + "proc-macro2", +] + +[[package]] +name = "rayon" +version = "1.10.0" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "b418a60154510ca1a002a752ca9714984e21e4241e804d32555251faf8b78ffa" +dependencies = [ + "either", + "rayon-core", +] + +[[package]] +name = "rayon-core" +version = "1.12.1" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "1465873a3dfdaa8ae7cb14b4383657caab0b3e8a0aa9ae8e04b044854c8dfce2" +dependencies = [ + "crossbeam-deque", + "crossbeam-utils", +] + +[[package]] +name = "redox_syscall" +version = "0.5.8" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "03a862b389f93e68874fbf580b9de08dd02facb9a788ebadaf4a3fd33cf58834" +dependencies = [ + "bitflags", +] + +[[package]] +name = "regex" +version = "1.11.1" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "b544ef1b4eac5dc2db33ea63606ae9ffcfac26c1416a2806ae0bf5f56b201191" +dependencies = [ + "aho-corasick", + "memchr", + "regex-automata", + "regex-syntax", +] + +[[package]] +name = "regex-automata" +version = "0.4.9" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "809e8dc61f6de73b46c85f4c96486310fe304c434cfa43669d7b40f711150908" +dependencies = [ + "aho-corasick", + "memchr", + "regex-syntax", +] + +[[package]] +name = "regex-syntax" +version = "0.8.5" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "2b15c43186be67a4fd63bee50d0303afffcef381492ebe2c5d87f324e1b8815c" + +[[package]] +name = "rustc-hash" +version = "1.1.0" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "08d43f7aa6b08d49f382cde6a7982047c3426db949b1424bc4b7ec9ae12c6ce2" + +[[package]] +name = "rustix" +version = "0.38.44" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "fdb5bc1ae2baa591800df16c9ca78619bf65c0488b41b96ccec5d11220d8c154" +dependencies = [ + "bitflags", + "errno", + "libc", + "linux-raw-sys", + "windows-sys", +] + +[[package]] +name = "scopeguard" +version = "1.2.0" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "94143f37725109f92c262ed2cf5e59bce7498c01bcc1502d7b9afe439a4e9f49" + +[[package]] +name = "shlex" +version = "1.3.0" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "0fda2ff0d084019ba4d7c6f371c95d8fd75ce3524c3cb8fb653a3023f6323e64" + +[[package]] +name = "smallvec" +version = "1.13.2" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "3c5e1a9a646d36c3599cd173a41282daf47c44583ad367b8e6837255952e5c67" + +[[package]] +name = "syn" +version = "2.0.98" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "36147f1a48ae0ec2b5b3bc5b537d267457555a10dc06f3dbc8cb11ba3006d3b1" +dependencies = [ + "proc-macro2", + "quote", + "unicode-ident", +] + +[[package]] +name = "target-lexicon" +version = "0.12.16" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "61c41af27dd6d1e27b1b16b489db798443478cef1f06a660c96db617ba5de3b1" + +[[package]] +name = "thiserror" +version = "2.0.12" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "567b8a2dae586314f7be2a752ec7474332959c6460e02bde30d702a66d488708" +dependencies = [ + "thiserror-impl", +] + +[[package]] +name = "thiserror-impl" +version = "2.0.12" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "7f7cf42b4507d8ea322120659672cf1b9dbb93f8f2d4ecfd6e51350ff5b17a1d" +dependencies = [ + "proc-macro2", + "quote", + "syn", +] + +[[package]] +name = "unicode-ident" +version = "1.0.16" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "a210d160f08b701c8721ba1c726c11662f877ea6b7094007e1ca9a1041945034" + +[[package]] +name = "unindent" +version = "0.2.3" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "c7de7d73e1754487cb58364ee906a499937a0dfabd86bcb980fa99ec8c8fa2ce" + +[[package]] +name = "which" +version = "4.4.2" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "87ba24419a2078cd2b0f2ede2691b6c66d8e47836da3b6db8265ebad47afbfc7" +dependencies = [ + "either", + "home", + "once_cell", + "rustix", +] + +[[package]] +name = "windows-sys" +version = "0.59.0" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "1e38bc4d79ed67fd075bcc251a1c39b32a1776bbe92e5bef1f0bf1f8c531853b" +dependencies = [ + "windows-targets", +] + +[[package]] +name = "windows-targets" +version = "0.52.6" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "9b724f72796e036ab90c1021d4780d4d3d648aca59e491e6b98e725b84e99973" +dependencies = [ + "windows_aarch64_gnullvm", + "windows_aarch64_msvc", + "windows_i686_gnu", + "windows_i686_gnullvm", + "windows_i686_msvc", + "windows_x86_64_gnu", + "windows_x86_64_gnullvm", + "windows_x86_64_msvc", +] + +[[package]] +name = "windows_aarch64_gnullvm" +version = "0.52.6" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "32a4622180e7a0ec044bb555404c800bc9fd9ec262ec147edd5989ccd0c02cd3" + +[[package]] +name = "windows_aarch64_msvc" +version = "0.52.6" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "09ec2a7bb152e2252b53fa7803150007879548bc709c039df7627cabbd05d469" + +[[package]] +name = "windows_i686_gnu" +version = "0.52.6" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "8e9b5ad5ab802e97eb8e295ac6720e509ee4c243f69d781394014ebfe8bbfa0b" + +[[package]] +name = "windows_i686_gnullvm" +version = "0.52.6" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "0eee52d38c090b3caa76c563b86c3a4bd71ef1a819287c19d586d7334ae8ed66" + +[[package]] +name = "windows_i686_msvc" +version = "0.52.6" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "240948bc05c5e7c6dabba28bf89d89ffce3e303022809e73deaefe4f6ec56c66" + +[[package]] +name = "windows_x86_64_gnu" +version = "0.52.6" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "147a5c80aabfbf0c7d901cb5895d1de30ef2907eb21fbbab29ca94c5b08b1a78" + +[[package]] +name = "windows_x86_64_gnullvm" +version = "0.52.6" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "24d5b23dc417412679681396f2b49f3de8c1473deb516bd34410872eff51ed0d" + +[[package]] +name = "windows_x86_64_msvc" +version = "0.52.6" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "589f6da84c646204747d1270a2a5661ea66ed1cced2631d546fdfb155959f9ec" diff --git a/polytope_feature/datacube/quadtree/Cargo.toml b/polytope_feature/datacube/quadtree/Cargo.toml new file mode 100644 index 000000000..8ac1653f0 --- /dev/null +++ b/polytope_feature/datacube/quadtree/Cargo.toml @@ -0,0 +1,16 @@ +[package] +name = "quadtree" +version = "0.1.0" +edition = "2021" + +[dependencies] +pyo3 = { version = "0.20", features = ["extension-module"] } +qhull = { path = "../../../../qhull-rs" } +rayon = "1.8" + +[lib] +crate-type = ["cdylib"] + +[profile.release] +debug = true + diff --git a/polytope_feature/datacube/quadtree/src/lib.rs b/polytope_feature/datacube/quadtree/src/lib.rs new file mode 100644 index 000000000..68f25101e --- /dev/null +++ b/polytope_feature/datacube/quadtree/src/lib.rs @@ -0,0 +1,736 @@ +#![allow(dead_code)] + +use pyo3::prelude::*; // Do not use * for importing here + +use qhull::Qh; +use std::collections::HashSet; +use std::error::Error; +use pyo3::exceptions::PyRuntimeError; +use qhull::*; + +use std::cmp::Ordering::Equal; + +// TODO: look at rust built in arena + + + + +// TODO: can we not replace this QuadPoint by just the index of the point in the list potentially? + + +#[derive(Debug)] +#[derive(Clone)] +#[pyclass] +struct QuadTreeNode { + points: Option>, + children: Vec, + + #[pyo3(get, set)] + center: (f64, f64), + size: (f64, f64), + depth: i32, +} + +impl QuadTreeNode { + fn sizeof(&self) -> usize { + let mut size = size_of::(); // Base struct size + // Dynamically allocated memory for points + if let Some(points) = &self.points { + // let points_size: usize = points.capacity() * size_of::(); + let points_size: usize = points.len() * size_of::(); + size += points_size; + } + // Dynamically allocated memory for children + let children_size = self.children.len() * size_of::(); + size += children_size; + size + } +} + + +#[derive(Debug)] +#[pyclass] +struct QuadTree { + nodes: Vec, +} + +#[pymethods] +impl QuadTree { + #[new] + fn new() -> Self { + QuadTree { + nodes: Vec::new(), + } + } + + fn sizeof(&self) -> usize { + let mut size = size_of::(); // Base struct size (Vec metadata) + // Memory allocated for `nodes` (Vec) + let nodes_size: usize = self.nodes.len() * size_of::(); + size += nodes_size; + // Sum sizes of each QuadTreeNode (including their allocated memory) + for (_i, node) in self.nodes.iter().enumerate() { + let node_size = node.sizeof(); + size += node_size; + } + size + } + + fn build_point_tree(&mut self, points: Vec<(f64, f64)>) { + self.create_node((0.0,0.0), (180.0, 90.0), 0); + points.iter().enumerate().for_each(|(index, _p)| { + self.insert(index, 0, &points); + }); + } + + + fn query_polygon(&mut self, quadtree_points: Vec<(f64, f64)>, node_idx: usize, mut polygon_points: Option>) -> PyResult> { + // Simulating a function that returns a Result + let mut results: HashSet = HashSet::new(); + + let processed_quadtree_points = self.process_points(quadtree_points); + + // let mut processed_polygon_points = self.convert_polygon(polygon_points).as_mut(); + // let mut processed_polygon_points = polygon_points.into_iter().map(|(x, y)| [x,y]).collect().as_mut(); + + // let mut processed_polygon_points: Vec<[f64; 2]> = polygon_points + // .into_iter() + // .map(|(x, y)| [x, y]) + // .collect(); + + let mut processed_polygon_points: Option> = polygon_points + .take() // Takes ownership, leaving None in the original variable + .map(|pts| pts.into_iter().map(|(x, y)| [x, y]).collect()); + + // Now, you can get a mutable reference + // let processed_polygon_points_ref: &mut Vec<[f64; 2]> = &mut processed_polygon_points; + + let query_result: Result<(), Box> = self._query_polygon(&processed_quadtree_points, node_idx, processed_polygon_points.as_mut(), &mut results); + + query_result.map_err(|e| PyErr::new::(e.to_string()))?; + + Ok(results) + } + + // /// Get the center of a node + fn get_center(&self, index: usize) -> PyResult<(f64, f64)> { + let nodes = &self.nodes; + nodes.get(index).map(|n| n.center).ok_or_else(|| { + pyo3::exceptions::PyIndexError::new_err("Invalid node index") + }) + } + + fn quadrant_rectangle_points(&self, node_idx: usize) -> PyResult> { + let (cx, cy) = self.get_center(node_idx)?; // Propagate error if get_center fails + let (sx, sy) = self.get_size(node_idx)?; // Propagate error if get_size fails + + Ok(vec![ + [cx - sx, cy - sy], + [cx - sx, cy + sy], + [cx + sx, cy - sy], + [cx + sx, cy + sy], + ]) + } + + fn find_nodes_in(&mut self, node_idx: usize) -> Vec { + let mut results = Vec::new(); + self.collect_points(&mut results, node_idx); + results + } + + fn get_children_idxs(&self, index: usize) -> Vec { + self.nodes.get(index).map_or_else(Vec::new, |node| node.children.to_vec()) + } + + fn get_point_idxs(&self, node_idx: usize) -> Vec { + self.nodes.get(node_idx) + .and_then(|n| n.points.as_ref()) // Get points if node exists + .map_or_else(Vec::new, |points| points.iter().map(|p| *p).collect()) + } + +} + + +impl QuadTree { + + const MAX: usize = 3; + const MAX_DEPTH: i32 = 20; + + fn process_points(&self, points: Vec<(f64, f64)>) -> Vec<[f64;2]> { + points.into_iter().map(|(x, y)| [x,y]).collect() + } + + // fn convert_polygon(&self, points: Option>) -> Option<&mut Vec<[f64; 2]>> { + // points.map(|pts| { + // let mut converted: Vec<[f64; 2]> = pts.iter().map(|&(x, y)| [x, y]).collect(); + // &mut converted + // }) + // } + + fn convert_polygon(&self, points: Option>) -> Option> { + points.map(|pts| pts.into_iter().map(|(x, y)| [x, y]).collect()) + } + + + fn create_node(&mut self, center: (f64, f64), size: (f64, f64), depth: i32) -> usize { + let nodes = &mut self.nodes; + let index = nodes.len(); + nodes.push(QuadTreeNode { + points: None, + children: Vec::new(), + center, + size, + depth, + }); + index + } + + fn get_depth(&self, index: usize) -> i32 { + let nodes = &self.nodes; + nodes.get(index).map(|n| n.depth).expect("Index exists in QuadTree arena") + } + + fn get_points_length(&self, index: usize) -> usize{ + let nodes = &self.nodes; + if let Some(n) = nodes.get(index) { + let point_count = n.points.as_ref().map_or(0, |v| v.len()); + point_count + } else { + panic!("Index exists in QuadTree arena"); + } + } + + fn get_size(&self, index: usize) -> PyResult<(f64, f64)> { + let nodes = &self.nodes; + nodes.get(index).map(|n| n.size).ok_or_else(|| { + pyo3::exceptions::PyIndexError::new_err("Invalid node index") + }) + } + + fn add_point_to_node(&mut self, index: usize, point_idx: usize) { + if let Some(n) = self.nodes.get_mut(index) { + // If the node already has points, check for duplicates + if let Some(ref mut points) = n.points { + if !points.iter().any(|pt| *pt == point_idx) { + points.push(point_idx); + } + } + // If there are no points yet, initialize the vector + else { + n.points = Some(vec![point_idx]); + } + } + } + + fn insert(&mut self, pt_index: usize, node_idx: usize, pts_ref: &Vec<(f64, f64)>) { + // Avoid allocating a new vector, check children directly + if self.nodes[node_idx].children.is_empty() { + self.add_point_to_node(node_idx, pt_index); + + // Avoid multiple calls to `self.get_points_length(node_idx)` + let points_len = self.get_points_length(node_idx); + let depth = self.get_depth(node_idx); + + if points_len > Self::MAX && depth < Self::MAX_DEPTH { + self.split(node_idx, pts_ref); + // TODO: here, can remove the points attribute of the node with node_idx + self.nodes[node_idx].points = None; + } + } else { + self.insert_into_children(pt_index, node_idx, pts_ref); + } + } + + + fn insert_into_children(&mut self, pt_index: usize, node_idx: usize, pts_ref: &Vec<(f64, f64)>) { + let (x,y) = pts_ref.get(pt_index).unwrap(); + let (cx, cy) = self.get_center(node_idx).unwrap(); + let child_idxs = self.get_children_idxs(node_idx); + + if *x <= cx { + if *y <= cy { + self.insert(pt_index, child_idxs[0], pts_ref); + } + if *y >= cy { + self.insert(pt_index, child_idxs[1], pts_ref); + } + } + if *x >= cx { + if *y <= cy { + self.insert(pt_index, child_idxs[2], pts_ref); + } + if *y >= cy { + self.insert(pt_index, child_idxs[3], pts_ref); + } + } + } + + fn add_child(&mut self, node_idx: usize, center: (f64, f64), size: (f64, f64), depth: i32) { + let child_idx = self.create_node( center, size, depth); + if let Some(n) = self.nodes.get_mut(node_idx) { + n.children.push(child_idx); + } + } + + fn split(&mut self, node_idx: usize, pts_ref: &Vec<(f64, f64)>) { + let (w, h) = self.get_size(node_idx).unwrap(); + let (x_center, y_center) = self.get_center(node_idx).unwrap(); + let node_depth = self.get_depth(node_idx); + + let (hx, hy) = (w * 0.5, h * 0.5); + + let new_centers = [ + (x_center - hx, y_center - hy), + (x_center - hx, y_center + hy), + (x_center + hx, y_center - hy), + (x_center + hx, y_center + hy), + ]; + + // Add children + for ¢er in &new_centers { + self.add_child(node_idx, center, (hx, hy), node_depth + 1); + } + + // Minimize locking scope + let points = self.nodes.get_mut(node_idx).and_then(|n| n.points.take()); + + // Process points outside the lock + if let Some(points) = points { + for node in points { + self.insert_into_children(node, node_idx, pts_ref); + } + } + } + + + fn collect_points(&mut self, results: &mut Vec, node_idx: usize) { + // Lock the nodes once and avoid locking multiple times + let nodes = &self.nodes; + + // Start by collecting the points of the current node + if let Some(n) = nodes.get(node_idx) { + if let Some(points) = &n.points { + results.extend(points.iter().map(|point| point)); // Use extend for more efficient pushing + } + } + + // Collect points from child nodes + let mut stack = Vec::new(); // Use a stack to avoid recursion overhead + stack.push(node_idx); // Start from the current node + + while let Some(current_node_idx) = stack.pop() { + let child_idxs = self.get_children_idxs(current_node_idx); + for child_idx in child_idxs { + stack.push(child_idx); // Add children to stack for later processing + } + + // Collect points of the child node + if let Some(n) = nodes.get(current_node_idx) { + if let Some(points) = &n.points { + results.extend(points.iter().map(|point| point)); // Efficiently extend the result + } + } + } + } + + fn get_node_items(&self, node_idx: usize) -> Vec { + let nodes = &self.nodes; + nodes.get(node_idx) + .and_then(|n| n.points.as_ref()) // Get `points` only if node exists + .map_or_else(Vec::new, |points| points.iter().map(|p| *p).collect()) + } + + // fn _query_polygon( + // &mut self, + // quadtree_points: &Vec<(f64, f64)>, + // node_idx: usize, + // polygon_points: Option<&mut Vec<(f64, f64)>>, + // results: &mut HashSet, + // ) -> Result<(), Box> { + // if let Some(points) = polygon_points { + // // Sort points only once + // points.sort_unstable_by(|a, b| a.0.partial_cmp(&b.0).unwrap()); + + // if *points == self.quadrant_rectangle_points(node_idx)? { + // results.extend(self.find_nodes_in(node_idx)); + // } else { + // let children_idxs = self.get_children_idxs(node_idx); + // if !children_idxs.is_empty() { + // let quadtree_center = self.get_center(node_idx)?; + + // let (left_polygon, right_polygon) = slice_in_two(Some(points), quadtree_center.0, 0)?; + // let (q1_polygon, q2_polygon) = slice_in_two(left_polygon.as_ref(), quadtree_center.1, 1)?; + // let (q3_polygon, q4_polygon) = slice_in_two(right_polygon.as_ref(), quadtree_center.1, 1)?; + + // if let Some(mut poly) = q1_polygon { + // self._query_polygon(quadtree_points, children_idxs[0], Some(poly.as_mut()), results)?; + // } + // if let Some(mut poly) = q2_polygon { + // self._query_polygon(quadtree_points, children_idxs[1], Some(poly.as_mut()), results)?; + // } + // if let Some(mut poly) = q3_polygon { + // self._query_polygon(quadtree_points, children_idxs[2], Some(poly.as_mut()), results)?; + // } + // if let Some(mut poly) = q4_polygon { + // self._query_polygon(quadtree_points, children_idxs[3], Some(poly.as_mut()), results)?; + // } + + // } else { + // let filtered_nodes: Vec = self + // .get_point_idxs(node_idx) + // .into_iter() + // .filter(|&node| is_contained_in(quadtree_points[node], &points)) + // .collect(); + // results.extend(filtered_nodes); + // } + // } + // } + // Ok(()) + // } + + fn _query_polygon( + &mut self, + quadtree_points: &Vec<[f64; 2]>, // Updated type + node_idx: usize, + polygon_points: Option<&mut Vec<[f64; 2]>>, // Updated type + results: &mut HashSet, + ) -> Result<(), Box> { + if let Some(points) = polygon_points { + // Sort points based on the first coordinate + points.sort_unstable_by(|a, b| a[0].partial_cmp(&b[0]).unwrap()); + + if *points == self.quadrant_rectangle_points(node_idx)? { + results.extend(self.find_nodes_in(node_idx)); + } else { + let children_idxs = self.get_children_idxs(node_idx); + if !children_idxs.is_empty() { + let quadtree_center = self.get_center(node_idx)?; + + let (left_polygon, right_polygon) = slice_in_two(Some(points), quadtree_center.0, 0)?; + let (q1_polygon, q2_polygon) = slice_in_two(left_polygon.as_ref(), quadtree_center.1, 1)?; + let (q3_polygon, q4_polygon) = slice_in_two(right_polygon.as_ref(), quadtree_center.1, 1)?; + + if let Some(mut poly) = q1_polygon { + self._query_polygon(quadtree_points, children_idxs[0], Some(poly.as_mut()), results)?; + } + if let Some(mut poly) = q2_polygon { + self._query_polygon(quadtree_points, children_idxs[1], Some(poly.as_mut()), results)?; + } + if let Some(mut poly) = q3_polygon { + self._query_polygon(quadtree_points, children_idxs[2], Some(poly.as_mut()), results)?; + } + if let Some(mut poly) = q4_polygon { + self._query_polygon(quadtree_points, children_idxs[3], Some(poly.as_mut()), results)?; + } + } else { + let filtered_nodes: Vec = self + .get_point_idxs(node_idx) + .into_iter() + .filter(|&node| is_contained_in(quadtree_points[node], &points)) + .collect(); + results.extend(filtered_nodes); + } + } + } + Ok(()) + } +} + + + +#[pymodule] +fn quadtree(_py: Python, m: &PyModule) -> PyResult<()> { + m.add_class::()?; + m.add_class::()?; + Ok(()) +} + +// fn is_contained_in(point: (f64, f64), polygon_points: &Vec<(f64, f64)>) -> bool { +// let (min_y_val, max_y_val) = _slice_2D_vertical_extents(&polygon_points, point.0); +// if min_y_val <= point.1 && point.1 <= max_y_val { +// return true +// } +// return false +// } + +fn is_contained_in(point: [f64; 2], polygon_points: &Vec<[f64; 2]>) -> bool { + let (min_y_val, max_y_val) = _slice_2D_vertical_extents(&polygon_points, point[0]); + point[1] >= min_y_val && point[1] <= max_y_val +} + + +// fn _slice_2D_vertical_extents(polygon_points: &Vec<(f64, f64)>, val: f64) -> (f64, f64){ +// let intersects = _find_intersects(&polygon_points, 0, val); +// intersects.into_iter().fold( +// (f64::INFINITY, f64::NEG_INFINITY), +// |(min, max), (_, y)| { +// (min.min(y), max.max(y)) // Only track the Y-values +// }, +// ) +// } + +// fn _slice_2D_vertical_extents(polygon_points: &Vec<[f64; 2]>, val: f64) -> (f64, f64) { +// let intersects = _find_intersects(polygon_points, 0, val); +// intersects.into_iter().fold( +// (f64::INFINITY, f64::NEG_INFINITY), // Initialize as [min, max] +// |(min, max), (_, y)| { +// (min.min(y), max.max(y)) // Only track the Y-values +// }, +// ) +// } + +fn _slice_2D_vertical_extents(polygon_points: &Vec<[f64; 2]>, val: f64) -> (f64, f64) { + // Get the intersection points with the vertical slice at `val` + let intersects = _find_intersects(polygon_points, 0, val); + + // Calculate the vertical extents (min and max Y-values) from the intersection points + intersects.into_iter().fold( + (f64::INFINITY, f64::NEG_INFINITY), // Initialize as [min, max] + |(min, max), intersect| { + // Destructure to access Y (intersect[1] is the Y-value) + let y = intersect[1]; + (min.min(y), max.max(y)) // Track the Y-values only + }, + ) +} + + +// RESTRICTED TO 2D POINTS FOR NOW +fn _find_intersects(polytope_points: &Vec<[f64; 2]>, slice_axis_idx: usize, value: f64) -> Vec<[f64; 2]>{ + let mut intersects: Vec<[f64; 2]> = vec![]; + let above_slice: Vec<[f64; 2]> = polytope_points + .iter() + .filter(|&point| { + let value_to_compare = if slice_axis_idx == 0 { point[0] } else { point[1] }; + value_to_compare >= value + }) + .copied() // Convert `&(f64, f64)` to `(f64, f64)` + .collect(); + + let below_slice: Vec<[f64; 2]> = polytope_points + .iter() + .filter(|&point| { + let value_to_compare = if slice_axis_idx == 0 { point[0] } else { point[1] }; + value_to_compare <= value + }) + .copied() // Convert `&(f64, f64)` to `(f64, f64)` + .collect(); + + + for a in &above_slice { + for b in &below_slice { + // Edge is incident with the slice plane, store b in intersects + if a[0] == b[0] && slice_axis_idx == 0 || a[1] == b[1] && slice_axis_idx == 1 { + intersects.push(*b); + continue; + } + let interp_coeff = (value - if slice_axis_idx == 0 { b[0] } else { b[1] }) + / (if slice_axis_idx == 0 { a[0] - b[0] } else { a[1] - b[1] }); + + let intersect = lerp(*a, *b, interp_coeff); + intersects.push(intersect); + } + } + intersects +} + + +fn lerp(a: [f64; 2], b: [f64; 2], t: f64) -> [f64; 2] { + [ + b[0] + t * (a[0] - b[0]), // Linear interpolation for x + b[1] + t * (a[1] - b[1]), // Linear interpolation for y + ] +} + + +fn polygon_extents(polytope_points: &Vec<[f64;2]>, slice_axis_idx: usize) -> (f64, f64){ + let (min_val, max_val) = polytope_points.into_iter().fold( + (f64::INFINITY, f64::NEG_INFINITY), // Start with extreme values + |(min, max), polytope_point| { + let value = if slice_axis_idx == 0 { polytope_point[0] } else { polytope_point[1] }; // Select the correct axis + (min.min(value), max.max(value)) // Update min and max + }, + ); + (min_val, max_val) +} + +// fn slice_in_two( +// polytope_points: Option<&Vec<(f64, f64)>>, +// value: f64, +// slice_axis_idx: usize, +// ) -> Result<(Option>, Option>), QhullError> { +// // Directly return if no points exist +// if let Some(polytope_points) = polytope_points { +// // Calculate the extents and intersections once +// let (x_lower, x_upper) = polygon_extents(&polytope_points, slice_axis_idx); +// let intersects = _find_intersects(&polytope_points, slice_axis_idx, value); + +// // If no intersections, directly handle the boundary cases +// if intersects.is_empty() { +// return Ok(if x_upper <= value { +// (Some(polytope_points.clone()), None) +// } else if value < x_lower { +// (None, Some(polytope_points.clone())) +// } else { +// (None, None) // Should never happen +// }); +// } + +// // Instead of partitioning into two vectors, we manually filter and extend +// let mut left_points = Vec::with_capacity(polytope_points.len()); +// let mut right_points = Vec::with_capacity(polytope_points.len()); + +// for &(x, y) in polytope_points { +// let value_to_compare = if slice_axis_idx == 0 { x } else { y }; +// if value_to_compare <= value { +// left_points.push((x, y)); +// } else { +// right_points.push((x, y)); +// } +// } + +// // Extend both left and right with intersection points +// left_points.extend(&intersects); +// right_points.extend(&intersects); + +// // Convert the points into polygons using find_qhull_points +// let left_polygon = find_qhull_points(&left_points)?; +// let right_polygon = find_qhull_points(&right_points)?; + +// return Ok((left_polygon, right_polygon)); +// } + +fn slice_in_two( + polytope_points: Option<&Vec<[f64; 2]>>, // Ensure points are in [f64; 2] format + value: f64, + slice_axis_idx: usize, +) -> Result<(Option>, Option>), QhullError> { + // Directly return if no points exist + if let Some(polytope_points) = polytope_points { + // Calculate the extents and intersections once + let (x_lower, x_upper) = polygon_extents(&polytope_points, slice_axis_idx); + let intersects: Vec<[f64; 2]> = _find_intersects(&polytope_points, slice_axis_idx, value); + + // If no intersections, directly handle the boundary cases + if intersects.is_empty() { + return Ok(if x_upper <= value { + (Some(polytope_points.clone()), None) + } else if value < x_lower { + (None, Some(polytope_points.clone())) + } else { + (None, None) // Should never happen + }); + } + + // Instead of partitioning into two vectors, we manually filter and extend + let mut left_points: Vec<[f64; 2]> = Vec::with_capacity(polytope_points.len()); + let mut right_points: Vec<[f64; 2]> = Vec::with_capacity(polytope_points.len()); + + for &point in polytope_points { + let value_to_compare = point[slice_axis_idx]; + if value_to_compare <= value { + left_points.push(point); + } else { + right_points.push(point); + } + } + + // Extend both left and right with intersection points + left_points.extend(intersects.iter().cloned()); + right_points.extend(intersects.iter().cloned()); + + // Convert the points into polygons using find_qhull_points + let left_polygon = find_qhull_points(&left_points)?; + let right_polygon = find_qhull_points(&right_points)?; + + return Ok((left_polygon, right_polygon)); + } + // Return None for both left and right if no polytope_points provided + Ok((None, None)) +} + + +fn change_points_for_qhull(points: &[(f64, f64)]) -> Vec<[f64; 2]> { + let mut result = Vec::with_capacity(points.len()); // Preallocate + for &(x, y) in points { + result.push([x, y]); + } + result +} + + + + +use std::fmt; + +#[derive(Debug)] +enum QhullError { + FlatError, + OtherError(String), +} + +impl fmt::Display for QhullError { + fn fmt(&self, f: &mut fmt::Formatter) -> fmt::Result { + match self { + QhullError::FlatError => write!(f, "QHull error: flat or invalid geometry"), + QhullError::OtherError(msg) => write!(f, "QHull error: {}", msg), + } + } +} + +impl Error for QhullError {} + + +fn find_qhull_points(points: &Vec<[f64; 2]>) -> Result>, QhullError> { + + // let converted = change_points_for_qhull(&points); + // let qh_result = Qh::builder() + // .compute(true) + // .build_from_iter(converted); + let qh_result = Qh::builder() + .compute(true) + .build_from_iter(points.iter().copied()); + // let qh_result = QhBuilder::default().build_managed(2, converted).unwrap(); + + match qh_result { + Ok(qh) => { + let num_vertices = qh.num_vertices(); // Get total number of vertices + let mut all_qhull_vertices_: HashSet = HashSet::with_capacity(num_vertices); + let mut all_qhull_vertices: Vec = Vec::with_capacity(num_vertices); + + // Process each simplex only once + for simplex in qh.simplices() { + for v in simplex.vertices().unwrap().iter() { + if let Some(index) = v.index(&qh) { + if all_qhull_vertices_.insert(index) { + all_qhull_vertices.push(index); + } + } + } + } + + // Allocate memory for final points + let mut actual_qhull_points: Vec<[f64; 2]> = Vec::with_capacity(all_qhull_vertices.len()); + + // Fetch actual points + for &idx in &all_qhull_vertices { + if let Some(&point) = points.get(idx) { + actual_qhull_points.push(point); + } + } + + Ok(Some(actual_qhull_points)) + } + Err(e) => { + let error_msg = e.to_string(); // Convert the error to a string + + if error_msg.contains("is flat") || error_msg.contains("less than"){ + Ok(None) + } else { + println!("QHull Error: {}", error_msg); + return Err(QhullError::OtherError("QHull Error".to_string())); + } + } + } +} + diff --git a/polytope_feature/datacube/quadtree/src/lib_old.rs b/polytope_feature/datacube/quadtree/src/lib_old.rs new file mode 100644 index 000000000..b5f4d6558 --- /dev/null +++ b/polytope_feature/datacube/quadtree/src/lib_old.rs @@ -0,0 +1,307 @@ +#![allow(dead_code)] + +use std::collections::HashMap; +use std::convert::TryFrom; +use std::time::Instant; +use std::mem; +use pyo3::prelude::*; // Do not use * for importing here +use pyo3::types::PyList; +use pyo3::exceptions::PyIndexError; +use pyo3::types::PyIterator; + + + +#[derive(Debug)] +#[derive(Clone)] +#[pyclass] +struct QuadTreeNode { + points: Option>, + children: Vec, + + #[pyo3(get, set)] + center: (f64, f64), + size: (f64, f64), + depth: i32, +} + +#[derive(Debug)] +#[derive(Clone)] +// #[pyclass] +struct QuadPoint { + // #[pyo3(get, set)] + item: (f64, f64), + // #[pyo3(get, set)] + index: i64, +} + + +// #[pymethods] +impl QuadPoint { + // #[new] + fn new(item: (f64, f64), index: i64) -> Self { + Self {item, index} + } +} + + +#[pymethods] +impl QuadTreeNode { + + const MAX: i32 = 3; + const MAX_DEPTH: i32 = 20; + + #[new] + fn new(x: f64, y: f64, size: Option<(f64, f64)>, depth: Option) -> Self { + Self { + points: Some(Vec::new()), + children: Vec::new(), + center: (x,y), + size: size.unwrap_or((180.0, 90.0)), + depth: depth.unwrap_or(0), + } + } + + fn build_point_tree(&mut self, points: Vec<(f64, f64)>) { + for (index, p) in points.iter().enumerate(){ + self.insert(p, index.try_into().unwrap()); + } + } + + fn quadrant_rectangle_points(&self) -> Vec<(f64, f64)> { + let (cx, cy) = self.center; + let (sx, sy) = self.size; + + let mut rect_points: Vec<(f64, f64)> = Vec::new(); + + rect_points.push((cx + sx, cy + sy)); + rect_points.push((cx + sx, cy - sy)); + rect_points.push((cx - sx, cy + sy)); + rect_points.push((cx - sx, cy - sy)); + + rect_points + } + + /// Return the length of the vector + fn __len__(&self) -> usize { + self.children.len() + } + + #[getter] + fn children(&self, py: Python) -> Py { + let py_list = PyList::empty(py); + for child in &self.children { + // py_list.append(child.clone()).unwrap(); + // let py_child = PyCell::new(py, child.clone()).unwrap(); + let py_child = unsafe {PyCell::new(py, child).unwrap().into();}; + py_list.append(py_child).unwrap(); + } + py_list.into() + } + + // #[getter] + // fn children(&self) -> Vec> { + // self.children.clone() // Cloning `Py` is cheap since it’s just a reference count increment + // } + + // #[getter] + // fn children(&self, py: Python) -> Py { + // let py_list = PyList::empty(py); + // for child in &self.children { + // // let py_child = Py::new(py, child).unwrap(); // Use Py::new to avoid cloning + // let py_child = PyCell::new(py, child).unwrap().into(); // Use PyCell to wrap + // py_list.append(py_child).unwrap(); + // } + // py_list.into() + // } + + // #[getter] + // fn children<'py>(&self, py: Python<'py>) -> Py { + // let py_list = PyList::empty(py); + // for child in &self.children { + // let py_child: Py = PyCell::new(py, child) + // .map(|cell| cell.into()) // Convert PyCell to a Python object + // .unwrap_or_else(|_| py.None()); // Handle errors safely + // py_list.append(py_child).unwrap(); + // } + // py_list.into() + // } + + + #[getter] + fn points(&self, py: Python) -> Py { + let py_list = PyList::empty(py); + if let Some(ref pts) = self.points { + for point in pts { + // let py_point = PyCell::new(py, point.clone()).unwrap(); + // py_list.append(py_point).unwrap(); + py_list.append(point.index); + } + } + py_list.into() + } + + fn find_nodes_in(&self) -> Vec { + let mut results = Vec::new(); + self.collect_points(&mut results); + results + } +} + +impl QuadTreeNode { + + fn get_node_items(&self) -> Vec<&(f64, f64)> { + let mut node_items: Vec<&(f64, f64)> = vec![]; + if let Some(points) = &self.points { + for point in points { + node_items.push(&point.item); + } + } + node_items + } + + fn insert(&mut self, item: &(f64, f64), index: i64) { + if self.children.is_empty() { + let mut node = QuadPoint::new(*item, index); + let mut node_items: Vec<&(f64, f64)> = vec![]; + if let Some(points) = self.points.as_mut() { + for point in &mut *points { + node_items.push(&point.item); + } + if !node_items.contains(&item) { + points.push(node); + } + if i32::try_from(points.len()).unwrap() > Self::MAX && self.depth < Self::MAX_DEPTH { + self.split(); + } + } + } + else { + self.insert_into_children(item, index); + } + } + + fn insert_into_children(&mut self, item: &(f64, f64), index: i64) { + let (x, y) = *item; + let (cx, cy) = self.center; + + if x <= cx { + if y <= cy { + self.children[0].insert(item, index); + } + if y >= cy { + self.children[1].insert(item, index); + } + } + if x >= cx { + if y <= cy { + self.children[2].insert(item, index); + } + if y >= cy { + self.children[3].insert(item, index); + } + } + } + + fn split(&mut self) { + + let (w, h) = self.size; + let hx: f64 = w/2.0; + let hy: f64 = h/2.0; + let (x_center, y_center) = self.center; + let new_centers: Vec<(f64, f64)> = vec![ + (x_center - hx, y_center - hy), + (x_center - hx, y_center + hy), + (x_center + hx, y_center - hy), + (x_center + hx, y_center + hy), + ]; + + self.children = new_centers.into_iter().map(|(x, y)| QuadTreeNode::new(x, y, Some((hx, hy)), Some(self.depth + 1))).collect(); + + if let Some(points) = self.points.take() { + for node in points { + self.insert_into_children(&node.item, node.index); + } + } + } + + + /// **Recursive helper function (not exposed to Python)** + fn collect_points(&self, results: &mut Vec) { + if let Some(points) = &self.points { + for point in points { + results.push(point.index); // Clone values instead of using references + } + } + + for child in &self.children { + child.collect_points(results); + } + } + + fn heap_size(&self) -> usize { + let mut total_size = 0; + + // Count memory used by `points` vector + if let Some(points) = &self.points { + total_size += mem::size_of_val(points) + (points.len() * mem::size_of::()); + } + + // Count memory used by `children` vector + total_size += mem::size_of_val(&self.children) + (self.children.len() * mem::size_of::()); + + // Recursively add children's heap size + for child in &self.children { + total_size += child.heap_size(); + } + + total_size + } + +} + +#[pymodule] +fn quadtree(_py: Python, m: &PyModule) -> PyResult<()> { + // m.add_class::()?; + m.add_class::()?; + Ok(()) +} + + +fn main() { + let mut quadtree: QuadTreeNode = QuadTreeNode::new(0.0, 0.0, Some((180.0, 90.0)), Some(0)); + println!("{:?}", quadtree.quadrant_rectangle_points()); + quadtree.insert(&(50.0, 50.0), 0); + if let Some(ref points) = quadtree.points { + println!("{}", points.len()); + println!("{:?}", points[0].item); + } + + let dx = 0.1; + let dy= 0.25; + + let mut grid = Vec::new(); + + let mut x = -180.0; + while x <= 180.0 { + let mut y = -90.0; + while y <= 90.0 { + grid.push((x, y)); + y += dy; + } + x += dx; + } + + println!("{}", grid.len()); + let start = Instant::now(); + quadtree.build_point_tree(grid); + let duration = start.elapsed(); + // println!("{:?}", quadtree); + println!("Time taken: {:?}", duration); + println!("Memory of the tree: {:?}", quadtree.heap_size() + mem::size_of_val(&quadtree)); +} + + + +// NOTE: everything/ every method that has to do with Polytope objects can be implemented in Python +// as a standalone method which takes in the Python objects and the QuadPoint or QuadTreeNode + diff --git a/polytope_feature/datacube/quadtree_additional_operations.py b/polytope_feature/datacube/quadtree_additional_operations.py new file mode 100644 index 000000000..f34c9fd65 --- /dev/null +++ b/polytope_feature/datacube/quadtree_additional_operations.py @@ -0,0 +1,112 @@ +from quadtree import QuadTree + +from ..engine.hullslicer import slice +from ..engine.slicing_tools import slice_in_two + + +def is_contained_in(point, polygon): + # implement method to check if the node point is inside the polygon + node_x, node_y = point + + sliced_vertical_polygon = slice(polygon, polygon._axes[0], node_x, 0) + if sliced_vertical_polygon: + lower_y, upper_y, idx = sliced_vertical_polygon.extents(polygon._axes[1]) + if lower_y <= node_y <= upper_y: + return True + return False + + +# def query_polygon(quadtree_points, quadtree: QuadTree, node_idx, polygon, results=None): +# # intersect quad tree with a 2D polygon +# if results is None: +# results = set() + +# # intersect the children with the polygon +# # TODO: here, we create None polygons... think about how to handle them +# if polygon is None: +# return results +# else: +# polygon_points = {tuple(point) for point in polygon.points} +# # TODO: are these the right points which we are comparing, ie the points on the polygon +# # and the points on the rectangle quadrant? +# if sorted(list(polygon_points)) == quadtree.quadrant_rectangle_points(node_idx): +# results.update(quadtree.find_nodes_in(node_idx)) +# else: +# children_idxs = quadtree.get_children_idxs(node_idx) +# if len(children_idxs) > 0: +# # first slice vertically +# quadtree_center = quadtree.get_center(node_idx) +# left_polygon, right_polygon = slice_in_two(polygon, quadtree_center[0], 0) + +# # then slice horizontally +# # ie need to slice the left and right polygons each in two to have the 4 quadrant polygons + +# q1_polygon, q2_polygon = slice_in_two(left_polygon, quadtree_center[1], 1) +# q3_polygon, q4_polygon = slice_in_two(right_polygon, quadtree_center[1], 1) + +# # now query these 4 polygons further down the quadtree +# query_polygon(quadtree_points, quadtree, children_idxs[0], q1_polygon, results) +# query_polygon(quadtree_points, quadtree, children_idxs[1], q2_polygon, results) +# query_polygon(quadtree_points, quadtree, children_idxs[2], q3_polygon, results) +# query_polygon(quadtree_points, quadtree, children_idxs[3], q4_polygon, results) + +# # TODO: try optimisation: take bbox of polygon and quickly remove the results that are not in bbox already +# else: +# # print(quadtree.get_point_idxs(node_idx)) +# # print(polygon.points) +# # print(len(children_idxs)) +# results.update( +# node for node in quadtree.get_point_idxs(node_idx) if is_contained_in(quadtree_points[node], polygon) +# ) +# return results + + +def query_polygon(quadtree_points, quadtree, node_idx, polygon): + results = set() + _query_polygon(quadtree_points, quadtree, node_idx, polygon, results) + return results + + +def _query_polygon(quadtree_points, quadtree: QuadTree, node_idx, polygon, results): + # intersect quad tree with a 2D polygon + # if results is None: + # results = set() + + # intersect the children with the polygon + # TODO: here, we create None polygons... think about how to handle them + if polygon is None: + return results + else: + polygon_points = {tuple(point) for point in polygon.points} + # TODO: are these the right points which we are comparing, ie the points on the polygon + # and the points on the rectangle quadrant? + if sorted(list(polygon_points)) == quadtree.quadrant_rectangle_points(node_idx): + results.update(quadtree.find_nodes_in(node_idx)) + else: + children_idxs = quadtree.get_children_idxs(node_idx) + if len(children_idxs) > 0: + # first slice vertically + quadtree_center = quadtree.get_center(node_idx) + left_polygon, right_polygon = slice_in_two(polygon, quadtree_center[0], 0) + + # then slice horizontally + # ie need to slice the left and right polygons each in two to have the 4 quadrant polygons + + q1_polygon, q2_polygon = slice_in_two(left_polygon, quadtree_center[1], 1) + q3_polygon, q4_polygon = slice_in_two(right_polygon, quadtree_center[1], 1) + + # now query these 4 polygons further down the quadtree + _query_polygon(quadtree_points, quadtree, children_idxs[0], q1_polygon, results) + _query_polygon(quadtree_points, quadtree, children_idxs[1], q2_polygon, results) + _query_polygon(quadtree_points, quadtree, children_idxs[2], q3_polygon, results) + _query_polygon(quadtree_points, quadtree, children_idxs[3], q4_polygon, results) + + # TODO: try optimisation: take bbox of polygon and quickly remove the results that are not in bbox already + else: + # print(quadtree.get_point_idxs(node_idx)) + # print(polygon.points) + # print(len(children_idxs)) + results.update( + node for node in quadtree.get_point_idxs(node_idx) if is_contained_in(quadtree_points[node], polygon) + ) + # return results diff --git a/polytope_feature/datacube/tensor_index_tree.py b/polytope_feature/datacube/tensor_index_tree.py index 604aa7441..ba1f841ef 100644 --- a/polytope_feature/datacube/tensor_index_tree.py +++ b/polytope_feature/datacube/tensor_index_tree.py @@ -109,6 +109,7 @@ def add_value(self, value): self.values = tuple(new_values) def create_child(self, axis, value, next_nodes): + # TODO: what if we remove the next nodes here? node = TensorIndexTree(axis, (value,)) existing_child = self.find_child(node) if not existing_child: @@ -116,6 +117,12 @@ def create_child(self, axis, value, next_nodes): return (node, next_nodes) return (existing_child, next_nodes) + def create_grid_child(self, axes, values): + # TODO: what if we remove the next nodes here? + node = GridNode(axes, values) + self.add_child(node) + return node + @property def parent(self): return self._parent @@ -138,12 +145,9 @@ def is_root(self): def find_child(self, node): index = self.children.bisect_left(node) - if index >= len(self.children): - return None - child = self.children[index] - if not child == node: - return None - return child + if index < len(self.children) and self.children[index] == node: + return self.children[index] + return None def add_node_layer_after(self, ax_name, vals): ax = IntDatacubeAxis() @@ -226,3 +230,82 @@ def get_ancestors(self): ancestors.append(current_node) current_node = current_node.parent return ancestors[::-1] + + +class GridNode(object): + def __init__(self, axes, values=tuple()): + # NOTE: the values here is a tuple so we can hash it + self.values = values + self.children = None + self._parent = None + self.result = [] + self.axes = axes + self.ancestors = [] + self.indexes = [] + self.hidden = False + + @property + def parent(self): + return self._parent + + @parent.setter + def set_parent(self, node): + if self.parent is not None: + self.parent.children.remove(self) + self._parent = node + self._parent.children.add(self) + + def __setitem__(self, key, value): + setattr(self, key, value) + + def __getitem__(self, key): + return getattr(self, key) + + def __delitem__(self, key): + return delattr(self, key) + + def __lt__(self, other): + return self.values < other.values + + def _collect_leaf_nodes(self, leaves): + if not self.children: + leaves.append(self) + self.ancestors.append(self) + else: + if len(self.children) == 0: + leaves.append(self) + self.ancestors.append(self) + for n in self.children: + for ancestor in self.ancestors: + n.ancestors.append(ancestor) + if self.axes != TensorIndexTree.root: + n.ancestors.append(self) + n._collect_leaf_nodes(leaves) + + def pprint(self, level=0): + if self.axes[0].name == "root": + logging.debug("\n") + logging.debug("\t" * level + "\u21b3" + str(self)) + if not self.children: + logging.debug("\t" * (level + 1) + "\u21b3" + str(self.result)) + else: + for child in self.children: + if not child.hidden: + child.pprint(level + 1) + if len(self.children) == 0: + logging.debug("\t" * (level + 1) + "\u21b3" + str(self.result)) + + def flatten(self): + path = DatacubePath() + path[self.axes[0].name] = (self.values[0],) + path[self.axes[1].name] = (self.values[1],) + # ancestors = self.get_ancestors() + # for ancestor in ancestors: + # path[ancestor.axis.name] = ancestor.values + return path + + def __repr__(self): + if self.axes != "root": + return f"{(self.axes[0].name, self.axes[1].name)}={self.values}" + else: + return f"{self.axes}" diff --git a/polytope_feature/datacube/transformations/datacube_mappers/datacube_mappers.py b/polytope_feature/datacube/transformations/datacube_mappers/datacube_mappers.py index e28566e69..23e891398 100644 --- a/polytope_feature/datacube/transformations/datacube_mappers/datacube_mappers.py +++ b/polytope_feature/datacube/transformations/datacube_mappers/datacube_mappers.py @@ -86,8 +86,8 @@ def find_second_idx(self, first_val, second_val): def unmap_first_val_to_start_line_idx(self, first_val): return self._final_transformation.unmap_first_val_to_start_line_idx(first_val) - def unmap(self, first_val, second_val): - return self._final_transformation.unmap(first_val, second_val) + def unmap(self, first_val, second_val, unmapped_idx=None): + return self._final_transformation.unmap(first_val, second_val, unmapped_idx) def find_modified_indexes(self, indexes, path, datacube, axis): if axis.name == self._mapped_axes()[0]: @@ -105,10 +105,16 @@ def unmap_path_key(self, key_value_path, leaf_path, unwanted_path, axis): unwanted_path[axis.name] = unwanted_val if axis.name == self._mapped_axes()[1]: first_val = unwanted_path[self._mapped_axes()[0]] - unmapped_idx = [] - for val in value: - unmapped_idx.append(self.unmap(first_val, (val,))) - # unmapped_idx = self.unmap(first_val, value) + unmapped_idx = leaf_path.get("index", None) + if unmapped_idx is not None and len(unmapped_idx) > 0: + # for val in value: + # unmapped_idx.append(self.unmap(first_val, (val,))) + unmapped_idx = list(unmapped_idx) + else: + unmapped_idx = [] + for val in value: + unmapped_idx.append(self.unmap(first_val, (val,))) + # unmapped_idx = self.unmap(first_val, value, unmapped_idx) leaf_path.pop(self._mapped_axes()[0], None) key_value_path.pop(axis.name) key_value_path[self.old_axis] = unmapped_idx @@ -136,5 +142,6 @@ def unmap_tree_node(self, node, unwanted_path): "regular": "RegularGridMapper", "reduced_ll": "ReducedLatLonMapper", "local_regular": "LocalRegularGridMapper", + "irregular": "IrregularGridMapper", "healpix_nested": "NestedHealpixGridMapper", } diff --git a/polytope_feature/datacube/transformations/datacube_mappers/mapper_types/healpix.py b/polytope_feature/datacube/transformations/datacube_mappers/mapper_types/healpix.py index 6fdf9d54a..556701312 100644 --- a/polytope_feature/datacube/transformations/datacube_mappers/mapper_types/healpix.py +++ b/polytope_feature/datacube/transformations/datacube_mappers/mapper_types/healpix.py @@ -133,7 +133,7 @@ def unmap_first_val_to_start_line_idx(self, first_val): else: return idx - def unmap(self, first_val, second_val): + def unmap(self, first_val, second_val, unmapped_idx=None): tol = 1e-8 first_value = [i for i in self._first_axis_vals if first_val[0] - tol <= i <= first_val[0] + tol][0] first_idx = self._first_axis_vals.index(first_value) diff --git a/polytope_feature/datacube/transformations/datacube_mappers/mapper_types/healpix_nested.py b/polytope_feature/datacube/transformations/datacube_mappers/mapper_types/healpix_nested.py index c71f59f5b..36a69492a 100644 --- a/polytope_feature/datacube/transformations/datacube_mappers/mapper_types/healpix_nested.py +++ b/polytope_feature/datacube/transformations/datacube_mappers/mapper_types/healpix_nested.py @@ -137,7 +137,7 @@ def unmap_first_val_to_start_line_idx(self, first_val): else: return idx - def unmap(self, first_val, second_val): + def unmap(self, first_val, second_val, unmapped_idx=None): tol = 1e-8 first_value = [i for i in self._first_axis_vals if first_val[0] - tol <= i <= first_val[0] + tol][0] first_idx = self._first_axis_vals.index(first_value) diff --git a/polytope_feature/datacube/transformations/datacube_mappers/mapper_types/irregular.py b/polytope_feature/datacube/transformations/datacube_mappers/mapper_types/irregular.py new file mode 100644 index 000000000..0dad07552 --- /dev/null +++ b/polytope_feature/datacube/transformations/datacube_mappers/mapper_types/irregular.py @@ -0,0 +1,23 @@ +from ..datacube_mappers import DatacubeMapper + + +class IrregularGridMapper(DatacubeMapper): + # def __init__(self, base_axis, mapped_axes, resolution, local_area=[]): + def __init__(self, base_axis, mapped_axes, resolution, md5_hash=None, local_area=[], axis_reversed=None): + self._mapped_axes = mapped_axes + self._base_axis = base_axis + self._resolution = resolution + self._axis_reversed = False + self.compressed_grid_axes = [self._mapped_axes[1]] + if md5_hash is not None: + self.md5_hash = md5_hash + else: + self.md5_hash = _md5_hash.get(resolution, None) + + def unmap(self, first_val, second_val, unmapped_idx=None): + # TODO: But to unmap for the irregular grid, need the request tree + # Suppose we get the idx value somehow from the tree, as an idx input + return unmapped_idx[0] + + +_md5_hash = {} diff --git a/polytope_feature/datacube/transformations/datacube_mappers/mapper_types/local_regular.py b/polytope_feature/datacube/transformations/datacube_mappers/mapper_types/local_regular.py index df611250e..fef0cb1f2 100644 --- a/polytope_feature/datacube/transformations/datacube_mappers/mapper_types/local_regular.py +++ b/polytope_feature/datacube/transformations/datacube_mappers/mapper_types/local_regular.py @@ -81,7 +81,7 @@ def unmap_first_val_to_start_line_idx(self, first_val): first_idx = self._first_axis_vals.index(first_val) return first_idx * self.second_resolution - def unmap(self, first_val, second_val): + def unmap(self, first_val, second_val, unmapped_idx=None): tol = 1e-8 first_val = [i for i in self._first_axis_vals if first_val[0] - tol <= i <= first_val[0] + tol][0] first_idx = self._first_axis_vals.index(first_val) diff --git a/polytope_feature/datacube/transformations/datacube_mappers/mapper_types/octahedral.py b/polytope_feature/datacube/transformations/datacube_mappers/mapper_types/octahedral.py index e03a45b67..415a5ba70 100644 --- a/polytope_feature/datacube/transformations/datacube_mappers/mapper_types/octahedral.py +++ b/polytope_feature/datacube/transformations/datacube_mappers/mapper_types/octahedral.py @@ -7883,7 +7883,7 @@ def find_second_axis_idx(self, first_val, second_val): second_idx = int(second_val[0] / second_axis_spacing) return (first_idx, second_idx) - def unmap(self, first_val, second_val): + def unmap(self, first_val, second_val, unmapped_idx=None): (first_idx, second_idx) = self.find_second_axis_idx(first_val, second_val) octahedral_index = self.axes_idx_to_octahedral_idx(first_idx, second_idx) return octahedral_index diff --git a/polytope_feature/datacube/transformations/datacube_mappers/mapper_types/reduced_ll.py b/polytope_feature/datacube/transformations/datacube_mappers/mapper_types/reduced_ll.py index 391e1bbe0..036060c09 100644 --- a/polytope_feature/datacube/transformations/datacube_mappers/mapper_types/reduced_ll.py +++ b/polytope_feature/datacube/transformations/datacube_mappers/mapper_types/reduced_ll.py @@ -1504,7 +1504,7 @@ def find_second_idx(self, first_val, second_val): second_idx = bisect.bisect_left(second_axis_vals, second_val - tol) return second_idx - def unmap(self, first_val, second_val): + def unmap(self, first_val, second_val, unmapped_idx=None): tol = 1e-8 first_value = [i for i in self._first_axis_vals if first_val[0] - tol <= i <= first_val[0] + tol][0] first_idx = self._first_axis_vals.index(first_value) diff --git a/polytope_feature/datacube/transformations/datacube_mappers/mapper_types/regular.py b/polytope_feature/datacube/transformations/datacube_mappers/mapper_types/regular.py index f4dfb6486..7bac1ea92 100644 --- a/polytope_feature/datacube/transformations/datacube_mappers/mapper_types/regular.py +++ b/polytope_feature/datacube/transformations/datacube_mappers/mapper_types/regular.py @@ -61,7 +61,7 @@ def unmap_first_val_to_start_line_idx(self, first_val): first_idx = self._first_axis_vals.index(first_val) return first_idx * 4 * self._resolution - def unmap(self, first_val, second_val): + def unmap(self, first_val, second_val, unmapped_idx=None): tol = 1e-8 first_val = [i for i in self._first_axis_vals if first_val[0] - tol <= i <= first_val[0] + tol][0] first_idx = self._first_axis_vals.index(first_val) diff --git a/polytope_feature/engine/engine.py b/polytope_feature/engine/engine.py index c714db0b2..d9674aae1 100644 --- a/polytope_feature/engine/engine.py +++ b/polytope_feature/engine/engine.py @@ -1,19 +1,43 @@ +from abc import abstractmethod from typing import List from ..datacube.backends.datacube import Datacube +from ..datacube.datacube_axis import UnsliceableDatacubeAxis from ..datacube.tensor_index_tree import TensorIndexTree from ..shapes import ConvexPolytope class Engine: - def __init__(self): - pass + def __init__(self, engine_options=None): + if engine_options is None: + engine_options = {} + self.engine_options = engine_options + self.ax_is_unsliceable = {} + + self.axis_values_between = {} + self.sliced_polytopes = {} + self.remapped_vals = {} + self.compressed_axes = [] def extract(self, datacube: Datacube, polytopes: List[ConvexPolytope]) -> TensorIndexTree: + # Delegate to the right slicer that the axes within the polytopes need to use pass + def check_slicer(self, ax): + # Return the slicer instance if ax is sliceable. + # If the ax is unsliceable, return None. + if isinstance(ax, UnsliceableDatacubeAxis): + return None + slicer_type = self.engine_options[ax.name] + slicer = self.generate_slicer(slicer_type) + return slicer + @staticmethod def default(): from .hullslicer import HullSlicer return HullSlicer() + + @abstractmethod + def _build_branch(self, ax, node, datacube, next_nodes, api): + pass diff --git a/polytope_feature/engine/hullslicer.py b/polytope_feature/engine/hullslicer.py index 9a99682a6..2d3bd0c22 100644 --- a/polytope_feature/engine/hullslicer.py +++ b/polytope_feature/engine/hullslicer.py @@ -1,41 +1,19 @@ import math from copy import copy from itertools import chain -from typing import List import scipy.spatial -from ..datacube.backends.datacube import Datacube -from ..datacube.datacube_axis import UnsliceableDatacubeAxis -from ..datacube.tensor_index_tree import TensorIndexTree from ..shapes import ConvexPolytope -from ..utility.combinatorics import group, tensor_product from ..utility.exceptions import UnsliceableShapeError from ..utility.geometry import lerp -from ..utility.list_tools import argmax, argmin, unique +from ..utility.list_tools import argmax, argmin from .engine import Engine class HullSlicer(Engine): def __init__(self): - self.ax_is_unsliceable = {} - self.axis_values_between = {} - self.has_value = {} - self.sliced_polytopes = {} - self.remapped_vals = {} - self.compressed_axes = [] - - def _unique_continuous_points(self, p: ConvexPolytope, datacube: Datacube): - for i, ax in enumerate(p._axes): - mapper = datacube.get_mapper(ax) - if self.ax_is_unsliceable.get(ax, None) is None: - self.ax_is_unsliceable[ax] = isinstance(mapper, UnsliceableDatacubeAxis) - if self.ax_is_unsliceable[ax]: - break - for j, val in enumerate(p.points): - p.points[j][i] = mapper.to_float(mapper.parse(p.points[j][i])) - # Remove duplicate points - unique(p.points) + super().__init__() def _build_unsliceable_child(self, polytope, ax, node, datacube, lowers, next_nodes, slice_axis_idx): if not polytope.is_flat: @@ -109,9 +87,9 @@ def remap_values(self, ax, value): self.remapped_vals[(value, ax.name)] = remapped_val return remapped_val - def _build_sliceable_child(self, polytope, ax, node, datacube, values, next_nodes, slice_axis_idx): + def _build_sliceable_child(self, polytope, ax, node, datacube, values, next_nodes, slice_axis_idx, api): for i, value in enumerate(values): - if i == 0 or ax.name not in self.compressed_axes: + if i == 0 or ax.name not in api.compressed_axes: fvalue = ax.to_float(value) new_polytope = slice(polytope, ax.name, fvalue, slice_axis_idx) remapped_val = self.remap_values(ax, value) @@ -125,8 +103,8 @@ def _build_sliceable_child(self, polytope, ax, node, datacube, values, next_node remapped_val = self.remap_values(ax, value) child.add_value(remapped_val) - def _build_branch(self, ax, node, datacube, next_nodes): - if ax.name not in self.compressed_axes: + def _build_branch(self, ax, node, datacube, next_nodes, api): + if ax.name not in api.compressed_axes: parent_node = node.parent right_unsliced_polytopes = [] for polytope in node["unsliced_polytopes"]: @@ -137,7 +115,7 @@ def _build_branch(self, ax, node, datacube, next_nodes): lower, upper, slice_axis_idx = polytope.extents(ax.name) # here, first check if the axis is an unsliceable axis and directly build node if it is # NOTE: we should have already created the ax_is_unsliceable cache before - if self.ax_is_unsliceable[ax.name]: + if api.ax_is_unsliceable[ax.name]: self._build_unsliceable_child(polytope, ax, node, datacube, [lower], next_nodes, slice_axis_idx) else: values = self.find_values_between(polytope, ax, node, datacube, lower, upper) @@ -148,7 +126,7 @@ def _build_branch(self, ax, node, datacube, next_nodes): # we have iterated all polytopes and we can now remove the node if we need to if len(values) == 0 and len(node.children) == 0: node.remove_branch() - self._build_sliceable_child(polytope, ax, node, datacube, values, next_nodes, slice_axis_idx) + self._build_sliceable_child(polytope, ax, node, datacube, values, next_nodes, slice_axis_idx, api) else: all_values = [] all_lowers = [] @@ -164,12 +142,12 @@ def _build_branch(self, ax, node, datacube, next_nodes): lower, upper, slice_axis_idx = polytope.extents(ax.name) if not first_slice_axis_idx: first_slice_axis_idx = slice_axis_idx - if self.ax_is_unsliceable[ax.name]: + if api.ax_is_unsliceable[ax.name]: all_lowers.append(lower) else: values = self.find_values_between(polytope, ax, node, datacube, lower, upper) all_values.extend(values) - if self.ax_is_unsliceable[ax.name]: + if api.ax_is_unsliceable[ax.name]: self._build_unsliceable_child( first_polytope, ax, node, datacube, all_lowers, next_nodes, first_slice_axis_idx ) @@ -177,75 +155,75 @@ def _build_branch(self, ax, node, datacube, next_nodes): if len(all_values) == 0: node.remove_branch() self._build_sliceable_child( - first_polytope, ax, node, datacube, all_values, next_nodes, first_slice_axis_idx + first_polytope, ax, node, datacube, all_values, next_nodes, first_slice_axis_idx, api ) del node["unsliced_polytopes"] - def find_compressed_axes(self, datacube, polytopes): - # First determine compressable axes from input polytopes - compressable_axes = [] - for polytope in polytopes: - if polytope.is_orthogonal: - for ax in polytope.axes(): - compressable_axes.append(ax) - # Cross check this list with list of compressable axis from datacube - # (should not include any merged or coupled axes) - for compressed_axis in compressable_axes: - if compressed_axis in datacube.compressed_axes: - self.compressed_axes.append(compressed_axis) - # add the last axis of the grid always (longitude) as a compressed axis - k, last_value = _, datacube.axes[k] = datacube.axes.popitem() - self.compressed_axes.append(k) - - def remove_compressed_axis_in_union(self, polytopes): - for p in polytopes: - if p.is_in_union: - for axis in p.axes(): - if axis == self.compressed_axes[-1]: - self.compressed_axes.remove(axis) - - def extract(self, datacube: Datacube, polytopes: List[ConvexPolytope]): - # Determine list of axes to compress - self.find_compressed_axes(datacube, polytopes) - - # remove compressed axes which are in a union - self.remove_compressed_axis_in_union(polytopes) - - # Convert the polytope points to float type to support triangulation and interpolation - for p in polytopes: - self._unique_continuous_points(p, datacube) - - groups, input_axes = group(polytopes) - datacube.validate(input_axes) - request = TensorIndexTree() - combinations = tensor_product(groups) - - # NOTE: could optimise here if we know combinations will always be for one request. - # Then we do not need to create a new index tree and merge it to request, but can just - # directly work on request and return it... - - for c in combinations: - r = TensorIndexTree() - new_c = [] - for combi in c: - if isinstance(combi, list): - new_c.extend(combi) - else: - new_c.append(combi) - r["unsliced_polytopes"] = set(new_c) - current_nodes = [r] - for ax in datacube.axes.values(): - next_nodes = [] - interm_next_nodes = [] - for node in current_nodes: - self._build_branch(ax, node, datacube, interm_next_nodes) - next_nodes.extend(interm_next_nodes) - interm_next_nodes = [] - current_nodes = next_nodes - - request.merge(r) - return request + # def find_compressed_axes(self, datacube, polytopes): + # # First determine compressable axes from input polytopes + # compressable_axes = [] + # for polytope in polytopes: + # if polytope.is_orthogonal: + # for ax in polytope.axes(): + # compressable_axes.append(ax) + # # Cross check this list with list of compressable axis from datacube + # # (should not include any merged or coupled axes) + # for compressed_axis in compressable_axes: + # if compressed_axis in datacube.compressed_axes: + # self.compressed_axes.append(compressed_axis) + # # add the last axis of the grid always (longitude) as a compressed axis + # k, last_value = _, datacube.axes[k] = datacube.axes.popitem() + # self.compressed_axes.append(k) + + # def remove_compressed_axis_in_union(self, polytopes): + # for p in polytopes: + # if p.is_in_union: + # for axis in p.axes(): + # if axis == self.compressed_axes[-1]: + # self.compressed_axes.remove(axis) + + # def extract(self, datacube: Datacube, polytopes: List[ConvexPolytope]): + # # Determine list of axes to compress + # self.find_compressed_axes(datacube, polytopes) + + # # remove compressed axes which are in a union + # self.remove_compressed_axis_in_union(polytopes) + + # # Convert the polytope points to float type to support triangulation and interpolation + # for p in polytopes: + # self._unique_continuous_points(p, datacube) + + # groups, input_axes = group(polytopes) + # datacube.validate(input_axes) + # request = TensorIndexTree() + # combinations = tensor_product(groups) + + # # NOTE: could optimise here if we know combinations will always be for one request. + # # Then we do not need to create a new index tree and merge it to request, but can just + # # directly work on request and return it... + + # for c in combinations: + # r = TensorIndexTree() + # new_c = [] + # for combi in c: + # if isinstance(combi, list): + # new_c.extend(combi) + # else: + # new_c.append(combi) + # r["unsliced_polytopes"] = set(new_c) + # current_nodes = [r] + # for ax in datacube.axes.values(): + # next_nodes = [] + # interm_next_nodes = [] + # for node in current_nodes: + # self._build_branch(ax, node, datacube, interm_next_nodes) + # next_nodes.extend(interm_next_nodes) + # interm_next_nodes = [] + # current_nodes = next_nodes + + # request.merge(r) + # return request def _find_intersects(polytope, slice_axis_idx, value): @@ -265,6 +243,9 @@ def _find_intersects(polytope, slice_axis_idx, value): # Linearly interpolate all coordinates of two points (a,b) of the polytope interp_coeff = (value - b[slice_axis_idx]) / (a[slice_axis_idx] - b[slice_axis_idx]) intersect = lerp(a, b, interp_coeff) + # print("HERE NOW LOOK") + # print(intersect) + # print(value) intersects.append(intersect) return intersects @@ -278,6 +259,7 @@ def _reduce_dimension(intersects, slice_axis_idx): def slice(polytope: ConvexPolytope, axis, value, slice_axis_idx): + # TODO: maybe these functions should go in the slicing tools? if polytope.is_flat: if value in chain(*polytope.points): intersects = [[value]] diff --git a/polytope_feature/engine/quadtree_slicer.py b/polytope_feature/engine/quadtree_slicer.py new file mode 100644 index 000000000..478bd2b01 --- /dev/null +++ b/polytope_feature/engine/quadtree_slicer.py @@ -0,0 +1,200 @@ +from copy import copy + +import numpy as np + +# from ..datacube.quad_tree import QuadTree +from quadtree import QuadTree + +from ..datacube.datacube_axis import IntDatacubeAxis +from ..datacube.quadtree_additional_operations import query_polygon +# from quadtree import query_polygon +from ..datacube.tensor_index_tree import TensorIndexTree +from .engine import Engine +import time + + +class QuadTreeSlicer(Engine): + def __init__(self, points): + # here need to construct quadtree, which is specific to datacube + # NOTE: should this be inside of the datacube instead that we create the quadtree? + # import time + + quad_tree = QuadTree() + print("START BUILDING QUAD TREE") + time0 = time.time() + points = [tuple(point) for point in points] + quad_tree.build_point_tree(points) + + print("SIZE OF THE QUAD TREE IS") + # print(quad_tree.sizeof()) + self.points = points + print("FINISH BUILDING QUAD TREE") + print(time.time() - time0) + self.quad_tree = quad_tree + + # method to slice polygon against quadtree + def extract(self, datacube, polytopes): + # need to find the points to extract within the polytopes (polygons here in 2D) + request = TensorIndexTree() + extracted_points = [] + for polytope in polytopes: + assert len(polytope._axes) == 2 + extracted_points.extend(self.extract_single(datacube, polytope)) + + # what data format do we return extracted points as? Append those points to the index tree? + + # NOTE: for now, we return the indices of the points in the point cloud, instead of lat/lon + for point in extracted_points: + # append each found leaf to the tree + idx = point + values_axis = IntDatacubeAxis() + values_axis.name = "values" + result = self.points[idx] + # TODO: make finding the axes objects nicer? + (child, _) = request.create_child(values_axis, idx, []) + child.result = result + + return request + + def extract_single(self, datacube, polytope): + # extract a single polygon + time1 = time.time() + # need to find points of the datacube contained within the polytope + # We do this by intersecting the datacube point cloud quad tree with the polytope here + # polygon_points = query_polygon(self.points, self.quad_tree, 0, polytope) + # polygon_points = query_polygon(self.points, self.quad_tree, 0, polytope) + polytope_points = [tuple(point) for point in polytope.points] + polygon_points = self.quad_tree.query_polygon(self.points, 0, polytope_points) + print("RUST QUERY POLYOGN TIME") + print(time.time() - time1) + return polygon_points + + def _build_branch(self, ax, node, datacube, next_nodes, api): + for polytope in node["unsliced_polytopes"]: + if ax.name in polytope._axes: + # here, first check if the axis is an unsliceable axis and directly build node if it is + + # NOTE: here, we only have sliceable children, since the unsliceable children are handled by the + # hullslicer engine? IS THIS TRUE? + self._build_sliceable_child(polytope, ax, node, datacube, next_nodes, api) + # TODO: what does this function actually return and what should it return? + # It just modifies the next_nodes? + del node["unsliced_polytopes"] + + def _build_sliceable_child(self, polytope, ax, node, datacube, next_nodes, api): + extracted_points = self.extract_single(datacube, polytope) + # TODO: add the sliced points as node to the tree and update the next_nodes + if len(extracted_points) == 0: + node.remove_branch() + + lat_ax = ax + lon_ax = datacube._axes["longitude"] + + for value in extracted_points: + # convert to float for slicing + lat_val = self.points[value][0] + lon_val = self.points[value][1] + + # store the native type + # (child, _) = node.create_child(lat_ax, lat_val, []) + # (grand_child, _) = child.create_child(lon_ax, lon_val, []) + + grand_child = node.create_grid_child((lat_ax, lon_ax), (lat_val, lon_val)) + # (grand_child, _) = child.create_child(lon_ax, lon_val, []) + # NOTE: the index of the point is stashed in the branches' result + grand_child.indexes = [value] + grand_child["unsliced_polytopes"] = copy(node["unsliced_polytopes"]) + grand_child["unsliced_polytopes"].remove(polytope) + # TODO: but now what happens to the second axis in the point cloud?? Do we create a second node for it?? + + +# from copy import copy + +# from ..datacube.datacube_axis import IntDatacubeAxis +# from ..datacube.quad_tree import QuadTree +# from ..datacube.tensor_index_tree import TensorIndexTree +# from .engine import Engine + + +# class QuadTreeSlicer(Engine): +# def __init__(self, points): +# # here need to construct quadtree, which is specific to datacube +# # NOTE: should this be inside of the datacube instead that we create the quadtree? +# import time +# quad_tree = QuadTree() +# print("START BUILDING QUAD TREE") +# time0 = time.time() +# quad_tree.build_point_tree(points) +# print("FINISH BUILDING QUAD TREE") +# print(time.time()-time0) +# self.quad_tree = quad_tree + +# # method to slice polygon against quadtree +# def extract(self, datacube, polytopes): +# # need to find the points to extract within the polytopes (polygons here in 2D) +# request = TensorIndexTree() +# extracted_points = [] +# for polytope in polytopes: +# assert len(polytope._axes) == 2 +# extracted_points.extend(self.extract_single(datacube, polytope)) + +# # what data format do we return extracted points as? Append those points to the index tree? + +# # NOTE: for now, we return the indices of the points in the point cloud, instead of lat/lon +# for point in extracted_points: +# # append each found leaf to the tree +# idx = point.index +# values_axis = IntDatacubeAxis() +# values_axis.name = "values" +# result = point.item +# # TODO: make finding the axes objects nicer? +# (child, _) = request.create_child(values_axis, idx, []) +# child.result = result + +# return request + +# def extract_single(self, datacube, polytope): +# # extract a single polygon + +# # need to find points of the datacube contained within the polytope +# # We do this by intersecting the datacube point cloud quad tree with the polytope here +# polygon_points = self.quad_tree.query_polygon(polytope) +# return polygon_points + +# def _build_branch(self, ax, node, datacube, next_nodes, api): +# for polytope in node["unsliced_polytopes"]: +# if ax.name in polytope._axes: +# # here, first check if the axis is an unsliceable axis and directly build node if it is + +# # NOTE: here, we only have sliceable children, since the unsliceable children are handled by the +# # hullslicer engine? IS THIS TRUE? +# self._build_sliceable_child(polytope, ax, node, datacube, next_nodes, api) +# # TODO: what does this function actually return and what should it return? +# # It just modifies the next_nodes? +# del node["unsliced_polytopes"] + +# def _build_sliceable_child(self, polytope, ax, node, datacube, next_nodes, api): +# extracted_points = self.extract_single(datacube, polytope) +# # TODO: add the sliced points as node to the tree and update the next_nodes +# if len(extracted_points) == 0: +# node.remove_branch() + +# for point in extracted_points: +# # convert to float for slicing +# value = point.index +# lat_val = point.item[0] +# lon_val = point.item[1] +# lat_ax = ax + +# # TODO: is there a nicer way to get this axis that does not depend on knowing +# # the axis name in advance? +# lon_ax = datacube._axes["longitude"] + +# # store the native type +# (child, _) = node.create_child(lat_ax, lat_val, []) +# (grand_child, _) = child.create_child(lon_ax, lon_val, []) +# # NOTE: the index of the point is stashed in the branches' result +# grand_child.indexes = [value] +# grand_child["unsliced_polytopes"] = copy(node["unsliced_polytopes"]) +# grand_child["unsliced_polytopes"].remove(polytope) +# # TODO: but now what happens to the second axis in the point cloud?? Do we create a second node for it?? diff --git a/polytope_feature/engine/slicing_tools.py b/polytope_feature/engine/slicing_tools.py new file mode 100644 index 000000000..14fa5312e --- /dev/null +++ b/polytope_feature/engine/slicing_tools.py @@ -0,0 +1,63 @@ +import scipy + +from ..shapes import ConvexPolytope +from .hullslicer import _find_intersects + + +def slice_in_two(polytope: ConvexPolytope, value, slice_axis_idx): + # print(value) + if polytope is None: + return (None, None) + else: + assert len(polytope.points[0]) == 2 + + x_lower, x_upper, _ = polytope.extents(polytope._axes[slice_axis_idx]) + + intersects = _find_intersects(polytope, slice_axis_idx, value) + + if len(intersects) == 0: + if x_upper <= value: + # The vertical slicing line does not intersect the polygon, which is on the left of the line + # So we keep the same polygon for now since it is unsliced + left_polygon = polytope + right_polygon = None + if value < x_lower: + left_polygon = None + right_polygon = polytope + else: + left_points = [p for p in polytope.points if p[slice_axis_idx] <= value] + right_points = [p for p in polytope.points if p[slice_axis_idx] >= value] + left_points.extend(intersects) + right_points.extend(intersects) + # find left polygon + try: + hull = scipy.spatial.ConvexHull(left_points) + vertices = hull.vertices + except scipy.spatial.qhull.QhullError as e: + if "less than" or "is flat" in str(e): + # NOTE: this happens when we slice a polygon that has a border which coincides with the quadrant + # line and we slice this additional border with the quadrant line again. + # This is not actually a polygon we want to consider so we ignore it + vertices = None + + if vertices is not None: + left_polygon = ConvexPolytope(polytope._axes, [left_points[i] for i in vertices]) + else: + left_polygon = None + + try: + hull = scipy.spatial.ConvexHull(right_points) + vertices = hull.vertices + except scipy.spatial.qhull.QhullError as e: + # NOTE: this happens when we slice a polygon that has a border which coincides with the quadrant + # line and we slice this additional border with the quadrant line again. + # This is not actually a polygon we want to consider so we ignore it + if "less than" or "is flat" in str(e): + vertices = None + + if vertices is not None: + right_polygon = ConvexPolytope(polytope._axes, [right_points[i] for i in vertices]) + else: + right_polygon = None + + return (left_polygon, right_polygon) diff --git a/polytope_feature/options.py b/polytope_feature/options.py index 2e88c9679..d309c0bf3 100644 --- a/polytope_feature/options.py +++ b/polytope_feature/options.py @@ -18,7 +18,7 @@ class CyclicConfig(TransformationConfig): class MapperConfig(TransformationConfig): name: Literal["mapper"] type: str = "" - resolution: Union[int, List[int]] = 0 + resolution: Optional[Union[int, List[int]]] = 0 axes: List[str] = [""] md5_hash: Optional[str] = None local: Optional[List[float]] = None @@ -61,7 +61,7 @@ class Config(ConfigModel): axis_config: List[AxisConfig] = [] compressed_axes_config: List[str] = [""] pre_path: Optional[Dict[str, path_subclasses_union]] = {} - alternative_axes: List[GribJumpAxesConfig] = [] + alternative_axes: Optional[List[GribJumpAxesConfig]] = [] class PolytopeOptions(ABC): diff --git a/polytope_feature/polytope.py b/polytope_feature/polytope.py index f18f10c4d..c787ceb00 100644 --- a/polytope_feature/polytope.py +++ b/polytope_feature/polytope.py @@ -1,9 +1,15 @@ -import logging from typing import List +from .datacube.backends.datacube import Datacube +from .datacube.datacube_axis import UnsliceableDatacubeAxis +from .datacube.tensor_index_tree import TensorIndexTree +from .engine.hullslicer import HullSlicer +from .engine.quadtree_slicer import QuadTreeSlicer from .options import PolytopeOptions from .shapes import ConvexPolytope +from .utility.combinatorics import group, tensor_product from .utility.exceptions import AxisOverdefinedError +from .utility.list_tools import unique class Request: @@ -39,33 +45,143 @@ def __repr__(self): class Polytope: - def __init__(self, datacube, engine=None, options=None, context=None): + def __init__( + self, + datacube, + options=None, + engine_options=None, + point_cloud_options=None, + context=None, + ): from .datacube import Datacube - from .engine import Engine if options is None: options = {} + if engine_options is None: + engine_options = {} - axis_options, compressed_axes_options, config, alternative_axes = PolytopeOptions.get_polytope_options(options) - + self.compressed_axes = [] self.context = context + axis_options, compressed_axes_options, config, alternative_axes = PolytopeOptions.get_polytope_options(options) self.datacube = Datacube.create( - datacube, config, axis_options, compressed_axes_options, alternative_axes, self.context + datacube, + config, + axis_options, + compressed_axes_options, + point_cloud_options, + alternative_axes, + self.context, ) - self.engine = engine if engine is not None else Engine.default() - self.time = 0 - - def slice(self, polytopes: List[ConvexPolytope]): + if engine_options == {}: + for ax_name in self.datacube._axes.keys(): + engine_options[ax_name] = "hullslicer" + self.engine_options = engine_options + self.engines = self.create_engines() + self.ax_is_unsliceable = {} + + def create_engines(self): + engines = {} + engine_types = set(self.engine_options.values()) + if "quadtree" in engine_types: + # TODO: need to get the corresponding point cloud from the datacube + quadtree_points = self.datacube.find_point_cloud() + engines["quadtree"] = QuadTreeSlicer(quadtree_points) + if "hullslicer" in engine_types: + engines["hullslicer"] = HullSlicer() + return engines + + def _unique_continuous_points(self, p: ConvexPolytope, datacube: Datacube): + for i, ax in enumerate(p._axes): + mapper = datacube.get_mapper(ax) + if self.ax_is_unsliceable.get(ax, None) is None: + self.ax_is_unsliceable[ax] = isinstance(mapper, UnsliceableDatacubeAxis) + if self.ax_is_unsliceable[ax]: + break + for j, val in enumerate(p.points): + p.points[j][i] = mapper.to_float(mapper.parse(p.points[j][i])) + # Remove duplicate points + unique(p.points) + + def slice(self, datacube, polytopes: List[ConvexPolytope]): """Low-level API which takes a polytope geometry object and uses it to slice the datacube""" - return self.engine.extract(self.datacube, polytopes) + self.find_compressed_axes(datacube, polytopes) - def retrieve(self, request: Request, method="standard"): + self.remove_compressed_axis_in_union(polytopes) + + # Convert the polytope points to float type to support triangulation and interpolation + for p in polytopes: + self._unique_continuous_points(p, datacube) + + groups, input_axes = group(polytopes) + datacube.validate(input_axes) + request = TensorIndexTree() + combinations = tensor_product(groups) + + # NOTE: could optimise here if we know combinations will always be for one request. + # Then we do not need to create a new index tree and merge it to request, but can just + # directly work on request and return it... + + for c in combinations: + r = TensorIndexTree() + new_c = [] + for combi in c: + if isinstance(combi, list): + new_c.extend(combi) + else: + new_c.append(combi) + r["unsliced_polytopes"] = set(new_c) + current_nodes = [r] + for ax in datacube.axes.values(): + engine = self.find_engine(ax) + next_nodes = [] + interm_next_nodes = [] + for node in current_nodes: + engine._build_branch(ax, node, datacube, interm_next_nodes, self) + next_nodes.extend(interm_next_nodes) + interm_next_nodes = [] + current_nodes = next_nodes + + request.merge(r) + return request + + def find_engine(self, ax): + slicer_type = self.engine_options[ax.name] + return self.engines[slicer_type] + + def old_retrieve(self, request: Request, method="standard"): """Higher-level API which takes a request and uses it to slice the datacube""" - logging.info("Starting request for %s ", self.context) - self.datacube.check_branching_axes(request) + # self.datacube.check_branching_axes(request) request_tree = self.engine.extract(self.datacube, request.polytopes()) - logging.info("Created request tree for %s ", self.context) - self.datacube.get(request_tree, self.context) - logging.info("Retrieved data for %s ", self.context) + self.datacube.get(request_tree) return request_tree + + def retrieve(self, request: Request, method="standard"): + """Higher-level API which takes a request and uses it to slice the datacube""" + request_tree = self.slice(self.datacube, request.polytopes()) + self.datacube.get(request_tree) + return request_tree + + def find_compressed_axes(self, datacube, polytopes): + # First determine compressable axes from input polytopes + compressable_axes = [] + for polytope in polytopes: + if polytope.is_orthogonal: + for ax in polytope.axes(): + compressable_axes.append(ax) + # Cross check this list with list of compressable axis from datacube + # (should not include any merged or coupled axes) + for compressed_axis in compressable_axes: + if compressed_axis in datacube.compressed_axes: + self.compressed_axes.append(compressed_axis) + + k, last_value = _, datacube.axes[k] = datacube.axes.popitem() + self.compressed_axes.append(k) + + def remove_compressed_axis_in_union(self, polytopes): + for p in polytopes: + if p.is_in_union: + for axis in p.axes(): + if axis in self.compressed_axes: + if axis == self.compressed_axes[-1]: + self.compressed_axes.remove(axis) diff --git a/polytope_feature/utility/engine_tools.py b/polytope_feature/utility/engine_tools.py new file mode 100644 index 000000000..08569ff4e --- /dev/null +++ b/polytope_feature/utility/engine_tools.py @@ -0,0 +1,31 @@ +from typing import List + +from ..datacube.backends.datacube import Datacube +from ..datacube.datacube_axis import UnsliceableDatacubeAxis +from ..datacube.tensor_index_tree import TensorIndexTree +from ..shapes import ConvexPolytope +from .combinatorics import group, tensor_product, unique + + +def unique_continuous_points_in_polytope(p: ConvexPolytope, datacube: Datacube): + # TODO: should this be in utility folder since it doesn't really depend on engine...? + for i, ax in enumerate(p._axes): + mapper = datacube.get_mapper(ax) + if isinstance(mapper, UnsliceableDatacubeAxis): + break + for j, val in enumerate(p.points): + p.points[j][i] = mapper.to_float(mapper.parse(p.points[j][i])) + # Remove duplicate points + unique(p.points) + + +def find_polytope_combinations(datacube: Datacube, polytopes: List[ConvexPolytope]) -> TensorIndexTree: + # here, we find the different possible polytope combinations that cover all of the datacube axes + + for p in polytopes: + unique_continuous_points_in_polytope(p, datacube) + + groups, input_axes = group(polytopes) + datacube.validate(input_axes) + combinations = tensor_product(groups) + return combinations diff --git a/profile.json.gz b/profile.json.gz new file mode 100644 index 000000000..34d84836c Binary files /dev/null and b/profile.json.gz differ diff --git a/pyproject.toml b/pyproject.toml index d9a2e1f4d..50a2dbbcb 100644 --- a/pyproject.toml +++ b/pyproject.toml @@ -5,4 +5,5 @@ requires = ["setuptools>=61.0"] build-backend = "setuptools.build_meta" [tool.pytest.ini_options] markers = ["internet: downloads test data from the internet (deselect with '-m \"not internet\"')", - "fdb: test which uses fdb (deselect with '-m \"not fdb\"')",] \ No newline at end of file + "fdb: test which uses fdb (deselect with '-m \"not fdb\"')",] +log_cli_level = "INFO" \ No newline at end of file diff --git a/tests/profiled_quadtree.profile b/tests/profiled_quadtree.profile new file mode 100644 index 000000000..597957624 Binary files /dev/null and b/tests/profiled_quadtree.profile differ diff --git a/tests/quadtree_slicer_profiler.py b/tests/quadtree_slicer_profiler.py new file mode 100644 index 000000000..ecb3f67a9 --- /dev/null +++ b/tests/quadtree_slicer_profiler.py @@ -0,0 +1,59 @@ +import cProfile + +import numpy as np +import pygribjump as gj + +from polytope_feature.polytope import Polytope, Request +from polytope_feature.shapes import Box + +options = { + "axis_config": [ + {"axis_name": "step", "transformations": [{"name": "type_change", "type": "int"}]}, + {"axis_name": "number", "transformations": [{"name": "type_change", "type": "int"}]}, + { + "axis_name": "date", + "transformations": [{"name": "merge", "other_axis": "time", "linkers": ["T", "00"]}], + }, + { + "axis_name": "values", + "transformations": [ + {"name": "mapper", "type": "irregular", "resolution": 1280, "axes": ["latitude", "longitude"]} + ], + }, + ], + "compressed_axes_config": [ + "longitude", + "latitude", + "levtype", + "step", + "date", + "domain", + "expver", + "param", + "class", + "stream", + "type", + ], + "pre_path": {"class": "od", "expver": "0001", "levtype": "sfc", "stream": "oper"}, +} +fdbdatacube = gj.GribJump() + +x = np.linspace(0, 100, 1000) +y = np.linspace(0, 100, 1000) +# create the mesh based on these arrays +X, Y = np.meshgrid(x, y) +X = X.reshape((np.prod(X.shape),)) +Y = Y.reshape((np.prod(Y.shape),)) +coords = zip(X, Y) +points = [list(coord) for coord in coords] +polytope = Box(["latitude", "longitude"], [1, 1], [20, 30]).polytope()[0] +API = Polytope( + request=Request(polytope), + datacube=fdbdatacube, + options=options, + engine_options={"latitude": "quadtree", "longitude": "quadtree"}, + point_cloud_options=points, +) +cProfile.runctx( + "API.engines['quadtree'].extract(API.datacube, [polytope])", globals(), locals(), "profiled_extract_quadtree.pstats" +) diff --git a/tests/test_bad_request_error.py b/tests/test_bad_request_error.py index b1154cf9c..dc2fc8fc2 100644 --- a/tests/test_bad_request_error.py +++ b/tests/test_bad_request_error.py @@ -60,6 +60,5 @@ def test_fdb_datacube(self): self.slicer = HullSlicer() self.API = Polytope( datacube=self.fdbdatacube, - engine=self.slicer, options=self.options, ) diff --git a/tests/test_cyclic_axis_over_negative_vals.py b/tests/test_cyclic_axis_over_negative_vals.py index 37bddd75f..5ba63c2e8 100644 --- a/tests/test_cyclic_axis_over_negative_vals.py +++ b/tests/test_cyclic_axis_over_negative_vals.py @@ -28,7 +28,7 @@ def setup_method(self, method): "compressed_axes_config": ["long", "level", "step", "date"], } self.slicer = HullSlicer() - self.API = Polytope(datacube=array, engine=self.slicer, options=options) + self.API = Polytope(datacube=array, options=options) # Testing different shapes diff --git a/tests/test_cyclic_axis_slicer_not_0.py b/tests/test_cyclic_axis_slicer_not_0.py index 84d36d9ce..b6ca580b8 100644 --- a/tests/test_cyclic_axis_slicer_not_0.py +++ b/tests/test_cyclic_axis_slicer_not_0.py @@ -2,7 +2,6 @@ import pandas as pd import xarray as xr -from polytope_feature.engine.hullslicer import HullSlicer from polytope_feature.polytope import Polytope, Request from polytope_feature.shapes import Box, Select @@ -27,10 +26,8 @@ def setup_method(self, method): ], "compressed_axes_config": ["long", "level", "step", "date"], } - self.slicer = HullSlicer() self.API = Polytope( datacube=array, - engine=self.slicer, options=self.options, ) diff --git a/tests/test_cyclic_axis_slicing.py b/tests/test_cyclic_axis_slicing.py index 3f8d53e8d..0f84c2c94 100644 --- a/tests/test_cyclic_axis_slicing.py +++ b/tests/test_cyclic_axis_slicing.py @@ -2,7 +2,6 @@ import pandas as pd import xarray as xr -from polytope_feature.engine.hullslicer import HullSlicer from polytope_feature.polytope import Polytope, Request from polytope_feature.shapes import Box, Select @@ -27,10 +26,8 @@ def setup_method(self, method): ], "compressed_axes_config": ["long", "level", "step", "date"], } - self.slicer = HullSlicer() self.API = Polytope( datacube=array, - engine=self.slicer, options=self.options, ) diff --git a/tests/test_cyclic_nearest.py b/tests/test_cyclic_nearest.py index 86e37c69c..7622eb915 100644 --- a/tests/test_cyclic_nearest.py +++ b/tests/test_cyclic_nearest.py @@ -3,7 +3,6 @@ from eccodes import codes_grib_find_nearest, codes_grib_new_from_file from helper_functions import download_test_data -from polytope_feature.engine.hullslicer import HullSlicer from polytope_feature.polytope import Polytope, Request from polytope_feature.shapes import Point, Select @@ -85,10 +84,8 @@ def test_regular_grid(self): Point(["latitude", "longitude"], [[39, 360 - 76.45]], method="nearest"), ) self.fdbdatacube = gj.GribJump() - self.slicer = HullSlicer() self.API = Polytope( datacube=self.fdbdatacube, - engine=self.slicer, options=self.options, ) result = self.API.retrieve(request) @@ -109,10 +106,8 @@ def test_regular_grid(self): Point(["latitude", "longitude"], [[39, -76.45]], method="nearest"), ) self.fdbdatacube = gj.GribJump() - self.slicer = HullSlicer() self.API = Polytope( datacube=self.fdbdatacube, - engine=self.slicer, options=self.options, ) result = self.API.retrieve(request) diff --git a/tests/test_cyclic_simple.py b/tests/test_cyclic_simple.py index d9974b7e7..538932a6b 100644 --- a/tests/test_cyclic_simple.py +++ b/tests/test_cyclic_simple.py @@ -2,7 +2,6 @@ import pandas as pd import xarray as xr -from polytope_feature.engine.hullslicer import HullSlicer from polytope_feature.polytope import Polytope, Request from polytope_feature.shapes import Box, Select @@ -27,10 +26,8 @@ def setup_method(self, method): ], "compressed_axes_config": ["long", "level", "step", "date"], } - self.slicer = HullSlicer() self.API = Polytope( datacube=array, - engine=self.slicer, options=options, ) diff --git a/tests/test_cyclic_snapping.py b/tests/test_cyclic_snapping.py index 599dc9573..5c55eaeea 100644 --- a/tests/test_cyclic_snapping.py +++ b/tests/test_cyclic_snapping.py @@ -21,7 +21,7 @@ def setup_method(self, method): "compressed_axes_config": ["long"], } self.slicer = HullSlicer() - self.API = Polytope(datacube=array, engine=self.slicer, options=options) + self.API = Polytope(datacube=array, options=options) # Testing different shapes diff --git a/tests/test_datacube_axes_init.py b/tests/test_datacube_axes_init.py index d92a08cf3..5dc6b329a 100644 --- a/tests/test_datacube_axes_init.py +++ b/tests/test_datacube_axes_init.py @@ -3,7 +3,6 @@ from helper_functions import download_test_data from polytope_feature.datacube.datacube_axis import FloatDatacubeAxis -from polytope_feature.engine.hullslicer import HullSlicer from polytope_feature.polytope import Polytope, Request from polytope_feature.shapes import Box, Select @@ -41,10 +40,8 @@ def setup_method(self, method): "type", ], } - self.slicer = HullSlicer() self.API = Polytope( datacube=latlon_array, - engine=self.slicer, options=self.options, ) self.datacube = self.API.datacube diff --git a/tests/test_datacube_xarray.py b/tests/test_datacube_xarray.py index f15a4f8bc..8f66eb822 100644 --- a/tests/test_datacube_xarray.py +++ b/tests/test_datacube_xarray.py @@ -24,7 +24,7 @@ def test_validate(self): array = xr.Dataset(data_vars=dict(param=(["x", "y", "z"], dims)), coords={"x": [1], "y": [1], "z": [1]}) array = array.to_array() - datacube = Datacube.create(array, axis_options={}) + datacube = Datacube.create(datacube=array, axis_options={}) datacube.validate(["x", "y", "z", "variable"]) datacube.validate(["x", "z", "y", "variable"]) @@ -53,7 +53,7 @@ def test_create(self): for d, v in array.coords.variables.items(): print(v.dtype) - datacube = Datacube.create(array, axis_options={}) + datacube = Datacube.create(datacube=array, axis_options={}) # Check the factory created the correct type of datacube assert isinstance(datacube, XArrayDatacube) diff --git a/tests/test_ecmwf_oper_data_fdb.py b/tests/test_ecmwf_oper_data_fdb.py index e848716c8..381cc566f 100644 --- a/tests/test_ecmwf_oper_data_fdb.py +++ b/tests/test_ecmwf_oper_data_fdb.py @@ -1,7 +1,6 @@ import pandas as pd import pytest -from polytope_feature.engine.hullslicer import HullSlicer from polytope_feature.polytope import Polytope, Request from polytope_feature.shapes import Box, Point, Select, Span @@ -59,10 +58,8 @@ def test_fdb_datacube(self): Box(["latitude", "longitude"], [0, 0], [0.2, 0.2]), ) self.fdbdatacube = gj.GribJump() - self.slicer = HullSlicer() self.API = Polytope( datacube=self.fdbdatacube, - engine=self.slicer, options=self.options, ) result = self.API.retrieve(request) @@ -87,10 +84,8 @@ def test_fdb_datacube_point(self): Point(["latitude", "longitude"], [[0.035149384216, 0.0]], method="surrounding"), ) self.fdbdatacube = gj.GribJump() - self.slicer = HullSlicer() self.API = Polytope( datacube=self.fdbdatacube, - engine=self.slicer, options=self.options, ) result = self.API.retrieve(request) @@ -116,10 +111,8 @@ def test_fdb_datacube_point_v2(self): Point(["latitude", "longitude"], [[0.035149384216, 0.0]], method="surrounding"), ) self.fdbdatacube = gj.GribJump() - self.slicer = HullSlicer() self.API = Polytope( datacube=self.fdbdatacube, - engine=self.slicer, options=self.options, ) result = self.API.retrieve(request) @@ -174,10 +167,8 @@ def test_fdb_datacube_point_step_not_compressed(self): Point(["latitude", "longitude"], [[0.035149384216, 0.0]], method="surrounding"), ) self.fdbdatacube = gj.GribJump() - self.slicer = HullSlicer() self.API = Polytope( datacube=self.fdbdatacube, - engine=self.slicer, options=self.options, ) result = self.API.retrieve(request) diff --git a/tests/test_engine_slicer.py b/tests/test_engine_slicer.py index 3392ccd3a..738b462eb 100644 --- a/tests/test_engine_slicer.py +++ b/tests/test_engine_slicer.py @@ -1,17 +1,17 @@ from polytope_feature.datacube.backends.mock import MockDatacube -from polytope_feature.engine.hullslicer import HullSlicer +from polytope_feature.polytope import Polytope from polytope_feature.shapes import Box, Polygon class TestEngineSlicer: def setup_method(self, method): - self.slicer = HullSlicer() + pass def test_2D_box(self): datacube = MockDatacube({"x": 100, "y": 100}) polytopes = Box(["x", "y"], lower_corner=[3, 3], upper_corner=[6, 6]).polytope() - result = self.slicer.extract(datacube, polytopes) - result.pprint() + API = Polytope(datacube=datacube, options={}) + result = API.slice(datacube, polytopes) assert len(result.leaves) == 4 total_leaves = 0 for leaf in result.leaves: @@ -21,8 +21,8 @@ def test_2D_box(self): def test_3D_box(self): datacube = MockDatacube({"x": 100, "y": 100, "z": 100}) polytopes = Box(["x", "y", "z"], lower_corner=[3, 3, 3], upper_corner=[6, 6, 6]).polytope() - result = self.slicer.extract(datacube, polytopes) - result.pprint() + API = Polytope(datacube=datacube, options={}) + result = API.slice(datacube, polytopes) assert len(result.leaves) == 4 * 4 total_leaves = 0 for leaf in result.leaves: @@ -32,8 +32,8 @@ def test_3D_box(self): def test_4D_box(self): datacube = MockDatacube({"x": 100, "y": 100, "z": 100, "q": 100}) polytopes = Box(["x", "y", "z", "q"], lower_corner=[3, 3, 3, 3], upper_corner=[6, 6, 6, 6]).polytope() - result = self.slicer.extract(datacube, polytopes) - result.pprint() + API = Polytope(datacube=datacube, options={}) + result = API.slice(datacube, polytopes) assert len(result.leaves) == 4 * 4 * 4 total_leaves = 0 for leaf in result.leaves: @@ -43,44 +43,40 @@ def test_4D_box(self): def test_triangle(self): datacube = MockDatacube({"x": 100, "y": 100}) triangle = Polygon(["x", "y"], [[3, 3], [3, 6], [6, 3]]).polytope() - result = self.slicer.extract(datacube, triangle) - result.pprint() - assert len(result.leaves) == 10 - # assert len(result.leaves) == 4 - # total_leaves = 0 - # for leaf in result.leaves: - # total_leaves += len(leaf.values) - # assert total_leaves == 4 + 3 + 2 + 1 + API = Polytope(datacube=datacube, options={}) + result = API.slice(datacube, triangle) + assert len(result.leaves) == 4 + 3 + 2 + 1 def test_reusable(self): datacube = MockDatacube({"x": 100, "y": 100}) polytopes = Polygon(["x", "y"], [[3, 3], [3, 6], [6, 3]]).polytope() - result = self.slicer.extract(datacube, polytopes) - result.pprint() - # assert len(result.leaves) == 4 - assert len(result.leaves) == 10 + API = Polytope(datacube=datacube, options={}) + result = API.slice(datacube, polytopes) + assert len(result.leaves) == 4 + 3 + 2 + 1 polytopes = Box(["x", "y"], lower_corner=[3, 3], upper_corner=[6, 6]).polytope() - result = self.slicer.extract(datacube, polytopes) - result.pprint() + result = API.slice(datacube, polytopes) assert len(result.leaves) == 4 def test_2D_box_get_function(self): datacube = MockDatacube({"x": 100, "y": 100}) polytopes = Box(["x", "y"], lower_corner=[2, -2], upper_corner=[4, -1]).polytope() - result = self.slicer.extract(datacube, polytopes) + API = Polytope(datacube=datacube, options={}) + result = API.slice(datacube, polytopes) datacube.get(result) result.pprint() def test_3D_box_get_function(self): datacube = MockDatacube({"x": 100, "y": 100, "z": 100}) polytopes = Box(["x", "y", "z"], lower_corner=[3, 2, -2], upper_corner=[6, 2, -1]).polytope() - result = self.slicer.extract(datacube, polytopes) + API = Polytope(datacube=datacube, options={}) + result = API.slice(datacube, polytopes) datacube.get(result) result.pprint() def test_3D_box_get_function2(self): datacube = MockDatacube({"x": 100, "y": 100, "z": 100}) polytopes = Box(["x", "y", "z"], lower_corner=[3, 2, 1], upper_corner=[6, 2, 1]).polytope() - result = self.slicer.extract(datacube, polytopes) + API = Polytope(datacube=datacube, options={}) + result = API.slice(datacube, polytopes) datacube.get(result) result.pprint() diff --git a/tests/test_fdb_datacube.py b/tests/test_fdb_datacube.py index 301ceec30..0b20ce78a 100644 --- a/tests/test_fdb_datacube.py +++ b/tests/test_fdb_datacube.py @@ -1,7 +1,6 @@ import pandas as pd import pytest -from polytope_feature.engine.hullslicer import HullSlicer from polytope_feature.polytope import Polytope, Request from polytope_feature.shapes import Box, Select, Span @@ -63,10 +62,8 @@ def test_fdb_datacube(self): Box(["latitude", "longitude"], [0, 0], [0.2, 0.2]), ) self.fdbdatacube = gj.GribJump() - self.slicer = HullSlicer() self.API = Polytope( datacube=self.fdbdatacube, - engine=self.slicer, options=self.options, ) result = self.API.retrieve(request) @@ -114,10 +111,8 @@ def test_fdb_datacube_select_grid(self): Span("longitude", 0, 0.070093457944), ) self.fdbdatacube = gj.GribJump() - self.slicer = HullSlicer() self.API = Polytope( datacube=self.fdbdatacube, - engine=self.slicer, options=self.options, ) result = self.API.retrieve(request) diff --git a/tests/test_fdb_unmap_tree.py b/tests/test_fdb_unmap_tree.py index a9dcf9cbe..c875433e7 100644 --- a/tests/test_fdb_unmap_tree.py +++ b/tests/test_fdb_unmap_tree.py @@ -1,7 +1,6 @@ import pandas as pd import pytest -from polytope_feature.engine.hullslicer import HullSlicer from polytope_feature.polytope import Polytope, Request from polytope_feature.shapes import Box, Select @@ -60,10 +59,8 @@ def test_fdb_datacube(self): Box(["latitude", "longitude"], [0, 0], [0.2, 0.2]), ) self.fdbdatacube = gj.GribJump() - self.slicer = HullSlicer() self.API = Polytope( datacube=self.fdbdatacube, - engine=self.slicer, options=self.options, ) result = self.API.retrieve(request) diff --git a/tests/test_float_type.py b/tests/test_float_type.py index 51df70df8..9916d6b1f 100644 --- a/tests/test_float_type.py +++ b/tests/test_float_type.py @@ -2,7 +2,6 @@ import pytest import xarray as xr -from polytope_feature.engine.hullslicer import HullSlicer from polytope_feature.polytope import Polytope, Request from polytope_feature.shapes import Select, Span @@ -19,9 +18,8 @@ def setup_method(self, method): "alt": np.arange(0.0, 20.0, 0.1), }, ) - self.slicer = HullSlicer() options = {"compressed_axes_config": ["lat", "long", "alt"]} - self.API = Polytope(datacube=array, engine=self.slicer, options=options) + self.API = Polytope(datacube=array, options=options) def test_slicing_span(self): request = Request(Span("lat", 4.1, 4.3), Select("long", [4.1]), Select("alt", [4.1])) diff --git a/tests/test_healpix_mapper.py b/tests/test_healpix_mapper.py index 8badbffc6..b2c9a7863 100644 --- a/tests/test_healpix_mapper.py +++ b/tests/test_healpix_mapper.py @@ -5,7 +5,6 @@ from polytope_feature.datacube.transformations.datacube_mappers.mapper_types.healpix import ( HealpixGridMapper, ) -from polytope_feature.engine.hullslicer import HullSlicer from polytope_feature.polytope import Polytope, Request from polytope_feature.shapes import Box, Select @@ -30,10 +29,8 @@ def setup_method(self, method): ], "compressed_axes_config": ["longitude", "latitude", "step", "time", "isobaricInhPa", "valid_time"], } - self.slicer = HullSlicer() self.API = Polytope( datacube=self.latlon_array, - engine=self.slicer, options=self.options, ) diff --git a/tests/test_hullslicer_engine.py b/tests/test_hullslicer_engine.py index 198b96b25..467ba3723 100644 --- a/tests/test_hullslicer_engine.py +++ b/tests/test_hullslicer_engine.py @@ -3,7 +3,6 @@ from polytope_feature.datacube.backends.xarray import XArrayDatacube from polytope_feature.datacube.tensor_index_tree import TensorIndexTree -from polytope_feature.engine.hullslicer import HullSlicer from polytope_feature.polytope import Polytope from polytope_feature.shapes import Box @@ -19,14 +18,13 @@ def setup_method(self, method): }, ) self.xarraydatacube = XArrayDatacube(self.array) - self.slicer = HullSlicer() options = {"compressed_axes_config": ["step", "level"]} - self.API = Polytope(datacube=self.array, engine=self.slicer, options=options) + self.API = Polytope(datacube=self.array, options=options) def test_extract(self): box = Box(["step", "level"], [3.0, 1.0], [6.0, 3.0]) polytope = box.polytope() - request = self.slicer.extract(self.xarraydatacube, polytope) + request = self.API.slice(self.xarraydatacube, polytope) assert request.axis == TensorIndexTree.root assert request.parent is None assert request.values is tuple() diff --git a/tests/test_icon_grid_unstructured.py b/tests/test_icon_grid_unstructured.py new file mode 100644 index 000000000..395e1b710 --- /dev/null +++ b/tests/test_icon_grid_unstructured.py @@ -0,0 +1,188 @@ +import math +import time + +import geopandas as gpd +import matplotlib.pyplot as plt +import numpy as np +import pandas as pd +import pytest +import xarray as xr +from earthkit import data +from helper_functions import find_nearest_latlon + +from polytope_feature.polytope import Polytope, Request +from polytope_feature.shapes import Box, Select + +# import iris + + + +class TestQuadTreeSlicer: + def setup_method(self, method): + self.engine_options = { + "step": "hullslicer", + "time": "hullslicer", + "heightAboveGround": "hullslicer", + "latitude": "quadtree", + "longitude": "quadtree", + } + print("SETTING UP THE XARRAY") + time_now = time.time() + + # ds = data.from_source( + # "file", "../../Downloads/icon-d2_germany_icosahedral_single-level_2025011000_024_2d_t_2m.grib2") + + # grid = xr.open_dataset("../../Downloads/icon_grid_0047_R19B07_L.nc", engine="netcdf4") + + ds = data.from_source("file", "../../Downloads/icon_global_icosahedral_single-level_2025011000_000_T_2M.grib2") + + grid = xr.open_dataset("../../Downloads/icon_grid_0026_R03B07_G.nc", engine="netcdf4") + + print(time.time() - time_now) + self.arr = ds.to_xarray(engine="cfgrib").t2m + + print(self.arr) + # print("AND GRID") + # print(grid.clon.values) + # print(len(grid.clon)) + self.longitudes = grid.clon.values * 180 / math.pi + self.latitudes = grid.clat.values * 180 / math.pi + print(grid) + + # self.arr = ds.votemper + + # self.latitudes = self.arr.latitude.values + # self.longitudes = self.arr.longitude.values + # self.latitudes = self.arr.coords["nav_lat"].values + # lats = [] + # for lat in self.latitudes: + # lats.extend(lat) + # self.latitudes = lats + # self.longitudes = self.arr.nav_lon.values + # lons = [] + # for lon in self.longitudes: + # lons.extend(lon) + # self.longitudes = lons + + # # self.arr["nav_lat_flat"] = (("grid_index",), self.arr.nav_lat.values.ravel()) + # # self.arr["nav_lon_flat"] = (("grid_index",), self.arr.nav_lon.values.ravel()) + + # # # Drop the x and y dimensions + # # self.arr = self.arr.drop_dims(["x", "y"]) + # nav_lat_flat = self.arr.nav_lat.values.ravel() + # nav_lon_flat = self.arr.nav_lon.values.ravel() + # deptht_flat = self.arr.deptht.values.ravel() + # interm_data = self.arr.data[0] + # new_interm_data = [] + # for data_pt in interm_data: + # new_data = data_pt.ravel() + # new_interm_data.append(new_data) + # print(len(interm_data)) + + # # Create a new dimension, for example, "grid_index" + # grid_index = np.arange(nav_lat_flat.size) + + # # Add the flattened `nav_lat` and `nav_lon` as variables + # # self.arr = self.arr.assign_coords(grid_index=("values", grid_index)) + # nav_lat_flat_da = xr.DataArray(nav_lat_flat, dims=["grid_index"], coords={"grid_index": grid_index}) + # nav_lon_flat_da = xr.DataArray(nav_lon_flat, dims=["grid_index"], coords={"grid_index": grid_index}) + + # # Drop x and y from the original DataArray + # # ds_cleaned = self.arr.drop(["x", "y"]) + + # # Combine everything into a new Dataset if needed + # # self.arr = xr.Dataset({ + # # "original_data": ds_cleaned, + # # "nav_lat_flat": nav_lat_flat_da, + # # "nav_lon_flat": nav_lon_flat_da, + # # }) + + # self.arr = xr.DataArray( + # new_interm_data, + # dims=("deptht", "values"), + # coords={ + # "deptht": deptht_flat, + # "values": grid_index, + # }, + # ) + + # print(self.arr) + # self.arr = self.arr.rename({"y": "lat", "nav_lon": "longitude", "x": "values"}) + self.points = list(zip(self.latitudes, self.longitudes)) + print((min(self.latitudes), max(self.latitudes), min(self.longitudes), max(self.longitudes))) + print("FINISH SETTING UP POINTS") + self.options = { + "axis_config": [ + { + "axis_name": "values", + "transformations": [ + {"name": "mapper", "type": "irregular", "resolution": 1280, "axes": ["latitude", "longitude"]} + ], + }, + ], + } + + @pytest.mark.fdb + def test_quad_tree_slicer_extract(self): + import datetime + + request = Request( + # Select("deptht", [0.5058], method="nearest"), + Select("time", [pd.Timestamp("2025-01-10")]), + Select("heightAboveGround", [2.0]), + Select("step", [datetime.timedelta(days=0)]), + # Select("time_counter", [pd.Timestamp("1979-02-15")]), + Box(["latitude", "longitude"], [0, 0], [10, 10]), + ) + + self.API = Polytope( + datacube=self.arr, + options=self.options, + engine_options=self.engine_options, + point_cloud_options=self.points, + ) + + time0 = time.time() + result = self.API.retrieve(request) + # result = self.API.slice(self.API.datacube, request.polytopes()) + time1 = time.time() + print("TIME TAKEN TO EXTRACT") + print(time1 - time0) + print(len(result.leaves)) + result.pprint() + + lats = [] + lons = [] + eccodes_lats = [] + eccodes_lons = [] + tol = 1e-8 + leaves = result.leaves + for i in range(len(leaves)): + cubepath = leaves[i].flatten() + lat = cubepath["latitude"][0] + lon = cubepath["longitude"][0] + lats.append(lat) + lons.append(lon) + + # # each variable in the netcdf file is a cube + # # cubes = iris.load('../../Downloads/votemper_ORAS5_1m_197902_grid_T_02.nc') + # # iris.save(cubes, '../../Downloads/votemper_ORAS5_1m_197902_grid_T_02.grib2') + # nearest_points = find_nearest_latlon( + # "../../Downloads/icon-d2_germany_icosahedral_single-level_2025011000_024_2d_t_2m.grib2", lat, lon) + # eccodes_lat = nearest_points[0][0]["lat"] + # eccodes_lon = nearest_points[0][0]["lon"] - 360 + # eccodes_lats.append(eccodes_lat) + # eccodes_lons.append(eccodes_lon) + # assert eccodes_lat - tol <= lat + # assert lat <= eccodes_lat + tol + # assert eccodes_lon - tol <= lon + # assert lon <= eccodes_lon + tol + + worldmap = gpd.read_file(gpd.datasets.get_path("naturalearth_lowres")) + fig, ax = plt.subplots(figsize=(12, 6)) + worldmap.plot(color="darkgrey", ax=ax) + + plt.scatter(lons, lats, s=18, c="red", cmap="YlOrRd") + plt.scatter(eccodes_lons, eccodes_lats, s=6, c="green") + plt.colorbar(label="Temperature") + plt.show() diff --git a/tests/test_icon_grid_unstructured_fdb.py b/tests/test_icon_grid_unstructured_fdb.py new file mode 100644 index 000000000..81762b891 --- /dev/null +++ b/tests/test_icon_grid_unstructured_fdb.py @@ -0,0 +1,152 @@ +import math + +# import iris +import os +import time + +import geopandas as gpd +import matplotlib.pyplot as plt +import numpy as np +import pandas as pd +import pytest +import xarray as xr +from earthkit import data +from helper_functions import find_nearest_latlon + +from polytope_feature.polytope import Polytope, Request +from polytope_feature.shapes import Box, Point, Select, Polygon + +os.environ["FDB_HOME"] = "/Users/male/git/fdb-new-home" + + +class TestQuadTreeSlicer: + def setup_method(self, method): + self.engine_options = { + "step": "hullslicer", + "date": "hullslicer", + "levtype": "hullslicer", + "param": "hullslicer", + "latitude": "quadtree", + "longitude": "quadtree", + } + print("SETTING UP THE XARRAY") + time_now = time.time() + + # ds = data.from_source( + # "file", "../../Downloads/icon-d2_germany_icosahedral_single-level_2025011000_024_2d_t_2m.grib2") + + # grid = xr.open_dataset("../../Downloads/icon_grid_0047_R19B07_L.nc", engine="netcdf4") + + ds = data.from_source("file", "../../Downloads/icon_global_icosahedral_single-level_2025011000_000_T_2M.grib2") + + grid = xr.open_dataset("../../Downloads/icon_grid_0026_R03B07_G.nc", engine="netcdf4") + + print(time.time() - time_now) + self.arr = ds.to_xarray(engine="cfgrib").t2m + + self.longitudes = grid.clon.values * 180 / math.pi + self.latitudes = grid.clat.values * 180 / math.pi + + self.points = list(zip(self.latitudes, self.longitudes)) + print((min(self.latitudes), max(self.latitudes), min(self.longitudes), max(self.longitudes))) + print("FINISH SETTING UP POINTS") + self.options = { + "axis_config": [ + {"axis_name": "step", "transformations": [{"name": "type_change", "type": "int"}]}, + { + "axis_name": "date", + "transformations": [{"name": "merge", "other_axis": "time", "linkers": ["T", "00"]}], + }, + { + "axis_name": "values", + "transformations": [ + { + "name": "mapper", + "type": "irregular", + "resolution": 0, + "axes": ["latitude", "longitude"], + "md5_hash": "f68071a8ac9bae4e965822afb963c04f", + } + ], + }, + ], + # "pre_path": {"time": "20250110", "heightAboveGround": "2"}, + "pre_path": {"date": "20250110"}, + } + + @pytest.mark.fdb + def test_quad_tree_slicer_extract(self): + import datetime + + import pygribjump as gj + + tri_side = 80 + triangle = Polygon(["latitude", "longitude"], [[0, tri_side], [0, 0], [tri_side, 0]]) + + request = Request( + # Select("deptht", [0.5058], method="nearest"), + Select("date", [pd.Timestamp("20250110T0000")]), + # Select("heightAboveGround", [2.0]), + # Select("step", [datetime.timedelta(days=0)]), + Select("step", [0]), + Select("param", ["167"]), + Select("levtype", ["sfc"]), + # Select("time_counter", [pd.Timestamp("1979-02-15")]), + # Box(["latitude", "longitude"], [0, 0], [5, 5]), + Box(["latitude", "longitude"], [0, 0], [80, 80]), + # triangle, + ) + + self.fdbdatacube = gj.GribJump() + + self.API = Polytope( + datacube=self.fdbdatacube, + options=self.options, + engine_options=self.engine_options, + point_cloud_options=self.points, + ) + + time0 = time.time() + result = self.API.retrieve(request) + # result = self.API.slice(self.API.datacube, request.polytopes()) + time1 = time.time() + print("TIME TAKEN TO EXTRACT") + print(time1 - time0) + print(len(result.leaves)) + result.pprint() + + lats = [] + lons = [] + eccodes_lats = [] + eccodes_lons = [] + tol = 1e-8 + leaves = result.leaves + for i in range(len(leaves)): + cubepath = leaves[i].flatten() + lat = cubepath["latitude"][0] + lon = cubepath["longitude"][0] + lats.append(lat) + lons.append(lon) + + # # each variable in the netcdf file is a cube + # # cubes = iris.load('../../Downloads/votemper_ORAS5_1m_197902_grid_T_02.nc') + # # iris.save(cubes, '../../Downloads/votemper_ORAS5_1m_197902_grid_T_02.grib2') + # nearest_points = find_nearest_latlon( + # "../../Downloads/icon-d2_germany_icosahedral_single-level_2025011000_024_2d_t_2m.grib2", lat, lon) + # eccodes_lat = nearest_points[0][0]["lat"] + # eccodes_lon = nearest_points[0][0]["lon"] - 360 + # eccodes_lats.append(eccodes_lat) + # eccodes_lons.append(eccodes_lon) + # assert eccodes_lat - tol <= lat + # assert lat <= eccodes_lat + tol + # assert eccodes_lon - tol <= lon + # assert lon <= eccodes_lon + tol + + # worldmap = gpd.read_file(gpd.datasets.get_path("naturalearth_lowres")) + fig, ax = plt.subplots(figsize=(12, 6)) + # worldmap.plot(color="darkgrey", ax=ax) + + plt.scatter(lons, lats, s=18, c="red", cmap="YlOrRd") + plt.scatter(eccodes_lons, eccodes_lats, s=6, c="green") + plt.colorbar(label="Temperature") + plt.show() diff --git a/tests/test_incomplete_tree_fdb.py b/tests/test_incomplete_tree_fdb.py index 2fbfcaf28..b388d37fd 100644 --- a/tests/test_incomplete_tree_fdb.py +++ b/tests/test_incomplete_tree_fdb.py @@ -2,7 +2,6 @@ import pytest from helper_functions import download_test_data -from polytope_feature.engine.hullslicer import HullSlicer from polytope_feature.polytope import Polytope, Request from polytope_feature.shapes import Select @@ -67,10 +66,8 @@ def test_incomplete_fdb_branch(self): Select("number", ["0"]), ) self.fdbdatacube = gj.GribJump() - self.slicer = HullSlicer() self.API = Polytope( datacube=self.fdbdatacube, - engine=self.slicer, options=self.options, ) result = self.API.retrieve(request) @@ -99,10 +96,8 @@ def test_incomplete_fdb_branch_2(self): Select("number", ["0"]), ) self.fdbdatacube = gj.GribJump() - self.slicer = HullSlicer() self.API = Polytope( datacube=self.fdbdatacube, - engine=self.slicer, options=self.options, ) result = self.API.retrieve(request) diff --git a/tests/test_local_grid_cyclic.py b/tests/test_local_grid_cyclic.py index 3f4bfd36c..32911dfce 100644 --- a/tests/test_local_grid_cyclic.py +++ b/tests/test_local_grid_cyclic.py @@ -1,7 +1,6 @@ import pandas as pd import pytest -from polytope_feature.engine.hullslicer import HullSlicer from polytope_feature.polytope import Polytope, Request from polytope_feature.shapes import Point, Select @@ -68,10 +67,8 @@ def test_fdb_datacube(self): Point(["latitude", "longitude"], [[-20, -20]]), ) self.fdbdatacube = gj.GribJump() - self.slicer = HullSlicer() self.API = Polytope( datacube=self.fdbdatacube, - engine=self.slicer, options=self.options, ) result = self.API.retrieve(request) @@ -97,10 +94,8 @@ def test_fdb_datacube_2(self): Point(["latitude", "longitude"], [[-20, 50 + 360]]), ) self.fdbdatacube = gj.GribJump() - self.slicer = HullSlicer() self.API = Polytope( datacube=self.fdbdatacube, - engine=self.slicer, options=self.options, ) result = self.API.retrieve(request) diff --git a/tests/test_local_regular_grid.py b/tests/test_local_regular_grid.py index 3f8f2d402..ac23bf172 100644 --- a/tests/test_local_regular_grid.py +++ b/tests/test_local_regular_grid.py @@ -1,7 +1,6 @@ import pandas as pd import pytest -from polytope_feature.engine.hullslicer import HullSlicer from polytope_feature.polytope import Polytope, Request from polytope_feature.shapes import Point, Select @@ -67,10 +66,8 @@ def test_fdb_datacube(self): Point(["latitude", "longitude"], [[0.16, 0.176]], method="nearest"), ) self.fdbdatacube = gj.GribJump() - self.slicer = HullSlicer() self.API = Polytope( datacube=self.fdbdatacube, - engine=self.slicer, options=self.options, ) result = self.API.retrieve(request) @@ -96,10 +93,8 @@ def test_point_outside_local_region(self): Point(["latitude", "longitude"], [[0.16, 61]], method="nearest"), ) self.fdbdatacube = gj.GribJump() - self.slicer = HullSlicer() self.API = Polytope( datacube=self.fdbdatacube, - engine=self.slicer, options=self.options, ) result = self.API.retrieve(request) @@ -125,10 +120,8 @@ def test_point_outside_local_region_2(self): Point(["latitude", "longitude"], [[41, 1]], method="nearest"), ) self.fdbdatacube = gj.GribJump() - self.slicer = HullSlicer() self.API = Polytope( datacube=self.fdbdatacube, - engine=self.slicer, options=self.options, ) result = self.API.retrieve(request) @@ -154,10 +147,8 @@ def test_point_outside_local_region_3(self): Point(["latitude", "longitude"], [[1, 61]]), ) self.fdbdatacube = gj.GribJump() - self.slicer = HullSlicer() self.API = Polytope( datacube=self.fdbdatacube, - engine=self.slicer, options=self.options, ) result = self.API.retrieve(request) @@ -182,10 +173,8 @@ def test_point_outside_local_region_4(self): Point(["latitude", "longitude"], [[41, 1]]), ) self.fdbdatacube = gj.GribJump() - self.slicer = HullSlicer() self.API = Polytope( datacube=self.fdbdatacube, - engine=self.slicer, options=self.options, ) result = self.API.retrieve(request) @@ -210,10 +199,8 @@ def test_point_outside_local_region_5(self): Point(["latitude", "longitude"], [[-41, 1]]), ) self.fdbdatacube = gj.GribJump() - self.slicer = HullSlicer() self.API = Polytope( datacube=self.fdbdatacube, - engine=self.slicer, options=self.options, ) result = self.API.retrieve(request) @@ -238,10 +225,8 @@ def test_point_outside_local_region_6(self): Point(["latitude", "longitude"], [[-30, -21]]), ) self.fdbdatacube = gj.GribJump() - self.slicer = HullSlicer() self.API = Polytope( datacube=self.fdbdatacube, - engine=self.slicer, options=self.options, ) result = self.API.retrieve(request) @@ -266,10 +251,8 @@ def test_point_outside_local_region_7(self): Point(["latitude", "longitude"], [[-41, 1]], method="nearest"), ) self.fdbdatacube = gj.GribJump() - self.slicer = HullSlicer() self.API = Polytope( datacube=self.fdbdatacube, - engine=self.slicer, options=self.options, ) result = self.API.retrieve(request) @@ -295,10 +278,8 @@ def test_point_outside_local_region_8(self): Point(["latitude", "longitude"], [[-30, -21]], method="nearest"), ) self.fdbdatacube = gj.GribJump() - self.slicer = HullSlicer() self.API = Polytope( datacube=self.fdbdatacube, - engine=self.slicer, options=self.options, ) result = self.API.retrieve(request) @@ -324,10 +305,8 @@ def test_point_outside_local_region_9(self): Point(["latitude", "longitude"], [[-30, -21]], method="surrounding"), ) self.fdbdatacube = gj.GribJump() - self.slicer = HullSlicer() self.API = Polytope( datacube=self.fdbdatacube, - engine=self.slicer, options=self.options, ) result = self.API.retrieve(request) diff --git a/tests/test_merge_octahedral_one_axis.py b/tests/test_merge_octahedral_one_axis.py index 4ddf90582..cde50d1d3 100644 --- a/tests/test_merge_octahedral_one_axis.py +++ b/tests/test_merge_octahedral_one_axis.py @@ -2,7 +2,6 @@ from earthkit import data from helper_functions import download_test_data -from polytope_feature.engine.hullslicer import HullSlicer from polytope_feature.polytope import Polytope, Request from polytope_feature.shapes import Box, Select @@ -28,10 +27,8 @@ def setup_method(self, method): ], "compressed_axes_config": ["longitude", "latitude", "surface", "step", "time", "valid_time", "number"], } - self.slicer = HullSlicer() self.API = Polytope( datacube=self.latlon_array, - engine=self.slicer, options=self.options, ) diff --git a/tests/test_merge_transformation.py b/tests/test_merge_transformation.py index 5d2311c35..99fc6d926 100644 --- a/tests/test_merge_transformation.py +++ b/tests/test_merge_transformation.py @@ -2,7 +2,6 @@ import pandas as pd import xarray as xr -from polytope_feature.engine.hullslicer import HullSlicer from polytope_feature.polytope import Polytope, Request from polytope_feature.shapes import Select @@ -27,10 +26,8 @@ def setup_method(self, method): ], "compressed_axes_config": ["date", "time"], } - self.slicer = HullSlicer() self.API = Polytope( datacube=self.array, - engine=self.slicer, options=self.options, ) diff --git a/tests/test_multiple_param_fdb.py b/tests/test_multiple_param_fdb.py index 11285af84..cd82d4731 100644 --- a/tests/test_multiple_param_fdb.py +++ b/tests/test_multiple_param_fdb.py @@ -1,7 +1,6 @@ import pandas as pd import pytest -from polytope_feature.engine.hullslicer import HullSlicer from polytope_feature.polytope import Polytope, Request from polytope_feature.shapes import Box, Select @@ -59,10 +58,8 @@ def test_fdb_datacube(self): Box(["latitude", "longitude"], [0, 0], [0.2, 0.2]), ) self.fdbdatacube = gj.GribJump() - self.slicer = HullSlicer() self.API = Polytope( datacube=self.fdbdatacube, - engine=self.slicer, options=self.options, ) result = self.API.retrieve(request) diff --git a/tests/test_octahedral_grid.py b/tests/test_octahedral_grid.py index f44ccf337..429794745 100644 --- a/tests/test_octahedral_grid.py +++ b/tests/test_octahedral_grid.py @@ -2,7 +2,6 @@ from earthkit import data from helper_functions import download_test_data, find_nearest_latlon -from polytope_feature.engine.hullslicer import HullSlicer from polytope_feature.polytope import Polytope, Request from polytope_feature.shapes import Box, Select @@ -27,10 +26,8 @@ def setup_method(self, method): ], "compressed_axes_config": ["longitude", "latitude", "number", "step", "time", "surface", "valid_time"], } - self.slicer = HullSlicer() self.API = Polytope( datacube=self.latlon_array, - engine=self.slicer, options=self.options, ) diff --git a/tests/test_orca_irregular_grid.py b/tests/test_orca_irregular_grid.py new file mode 100644 index 000000000..145af2079 --- /dev/null +++ b/tests/test_orca_irregular_grid.py @@ -0,0 +1,94 @@ +# import geopandas as gpd +# import matplotlib.pyplot as plt +import numpy as np +import pandas as pd +import pytest +from earthkit import data +from helper_functions import find_nearest_latlon + +from polytope_feature.polytope import Polytope, Request +from polytope_feature.shapes import Box, Select + + +class TestQuadTreeSlicer: + def setup_method(self, method): + self.engine_options = { + "step": "hullslicer", + "time": "hullslicer", + "latitude": "quadtree", + "longitude": "quadtree", + "oceanModelLayer": "hullslicer", + "valid_time": "hullslicer", + } + + ds = data.from_source("file", "../../Downloads/Reference_eORCA12_U_to_HEALPix_32.grib") + self.arr = ds.to_xarray(engine="cfgrib").avg_uox + + self.latitudes = self.arr.latitude.values + self.longitudes = self.arr.longitude.values + self.points = list(zip(self.latitudes, self.longitudes)) + self.options = { + "axis_config": [ + { + "axis_name": "values", + "transformations": [ + {"name": "mapper", "type": "irregular", "resolution": 1280, "axes": ["latitude", "longitude"]} + ], + }, + ], + } + + @pytest.mark.fdb + def test_quad_tree_slicer_extract(self): + request = Request( + Select("step", [np.timedelta64(0, "ns")]), + Select("oceanModelLayer", [65.0]), + Select("time", [pd.Timestamp("2017-09-06T00:00:00.000000000")]), + Box(["latitude", "longitude"], [65, 270], [75, 300]), + ) + + self.API = Polytope( + datacube=self.arr, + options=self.options, + engine_options=self.engine_options, + point_cloud_options=self.points, + ) + import time + + time0 = time.time() + result = self.API.retrieve(request) + time1 = time.time() + print("TIME TAKEN TO EXTRACT") + print(time1 - time0) + print(len(result.leaves)) + result.pprint() + + lats = [] + lons = [] + eccodes_lats = [] + eccodes_lons = [] + tol = 1e-8 + for i in range(len(result.leaves)): + cubepath = result.leaves[i].flatten() + lat = cubepath["latitude"][0] + lon = cubepath["longitude"][0] - 360 + lats.append(lat) + lons.append(lon) + nearest_points = find_nearest_latlon("../../Downloads/Reference_eORCA12_U_to_HEALPix_32.grib", lat, lon) + eccodes_lat = nearest_points[0][0]["lat"] + eccodes_lon = nearest_points[0][0]["lon"] - 360 + eccodes_lats.append(eccodes_lat) + eccodes_lons.append(eccodes_lon) + assert eccodes_lat - tol <= lat + assert lat <= eccodes_lat + tol + assert eccodes_lon - tol <= lon + assert lon <= eccodes_lon + tol + + # worldmap = gpd.read_file(gpd.datasets.get_path("naturalearth_lowres")) + # fig, ax = plt.subplots(figsize=(12, 6)) + # worldmap.plot(color="darkgrey", ax=ax) + + # plt.scatter(lons, lats, s=18, c="red", cmap="YlOrRd") + # plt.scatter(eccodes_lons, eccodes_lats, s=6, c="green") + # plt.colorbar(label="Temperature") + # plt.show() diff --git a/tests/test_orca_irregular_grid_seas5.py b/tests/test_orca_irregular_grid_seas5.py new file mode 100644 index 000000000..c8bc22fe5 --- /dev/null +++ b/tests/test_orca_irregular_grid_seas5.py @@ -0,0 +1,164 @@ +import time + +import geopandas as gpd +import matplotlib.pyplot as plt +import numpy as np +import pandas as pd +import pytest +import xarray as xr +from earthkit import data +from helper_functions import find_nearest_latlon + +from polytope_feature.polytope import Polytope, Request +from polytope_feature.shapes import Box, Select + +# import iris + + + +class TestQuadTreeSlicer: + def setup_method(self, method): + self.engine_options = { + "deptht": "hullslicer", + "latitude": "quadtree", + "longitude": "quadtree", + } + print("SETTING UP THE XARRAY") + time_now = time.time() + + # ds = data.from_source("file", "../../Downloads/votemper_ORAS5_1m_197902_grid_T_02.nc") + + ds = xr.open_dataset("../../Downloads/votemper_ORAS5_1m_197902_grid_T_02.nc", engine="netcdf4") + + print(time.time() - time_now) + # self.arr = ds.to_xarray(engine="cfgrib").avg_uox + self.arr = ds.votemper + + # self.latitudes = self.arr.latitude.values + # self.longitudes = self.arr.longitude.values + self.latitudes = self.arr.coords["nav_lat"].values + lats = [] + for lat in self.latitudes: + lats.extend(lat) + self.latitudes = lats + self.longitudes = self.arr.nav_lon.values + lons = [] + for lon in self.longitudes: + lons.extend(lon) + self.longitudes = lons + + # self.arr["nav_lat_flat"] = (("grid_index",), self.arr.nav_lat.values.ravel()) + # self.arr["nav_lon_flat"] = (("grid_index",), self.arr.nav_lon.values.ravel()) + + # # Drop the x and y dimensions + # self.arr = self.arr.drop_dims(["x", "y"]) + nav_lat_flat = self.arr.nav_lat.values.ravel() + nav_lon_flat = self.arr.nav_lon.values.ravel() + deptht_flat = self.arr.deptht.values.ravel() + interm_data = self.arr.data[0] + new_interm_data = [] + for data in interm_data: + new_data = data.ravel() + new_interm_data.append(new_data) + print(len(interm_data)) + + # Create a new dimension, for example, "grid_index" + grid_index = np.arange(nav_lat_flat.size) + + # Add the flattened `nav_lat` and `nav_lon` as variables + # self.arr = self.arr.assign_coords(grid_index=("values", grid_index)) + nav_lat_flat_da = xr.DataArray(nav_lat_flat, dims=["grid_index"], coords={"grid_index": grid_index}) + nav_lon_flat_da = xr.DataArray(nav_lon_flat, dims=["grid_index"], coords={"grid_index": grid_index}) + + # Drop x and y from the original DataArray + # ds_cleaned = self.arr.drop(["x", "y"]) + + # Combine everything into a new Dataset if needed + # self.arr = xr.Dataset({ + # "original_data": ds_cleaned, + # "nav_lat_flat": nav_lat_flat_da, + # "nav_lon_flat": nav_lon_flat_da, + # }) + + self.arr = xr.DataArray( + new_interm_data, + dims=("deptht", "values"), + coords={ + "deptht": deptht_flat, + "values": grid_index, + }, + ) + + print(self.arr) + # self.arr = self.arr.rename({"y": "lat", "nav_lon": "longitude", "x": "values"}) + self.points = list(zip(self.latitudes, self.longitudes)) + print("FINISH SETTING UP POINTS") + self.options = { + "axis_config": [ + { + "axis_name": "values", + "transformations": [ + {"name": "mapper", "type": "irregular", "resolution": 1280, "axes": ["latitude", "longitude"]} + ], + }, + ], + } + + @pytest.mark.fdb + def test_quad_tree_slicer_extract(self): + request = Request( + Select("deptht", [0.5058], method="nearest"), + # Select("time_counter", [pd.Timestamp("1979-02-15")]), + Box(["latitude", "longitude"], [-80, 50], [-50, 80]), + ) + + self.API = Polytope( + datacube=self.arr, + options=self.options, + engine_options=self.engine_options, + point_cloud_options=self.points, + ) + + time0 = time.time() + result = self.API.retrieve(request) + # result = self.API.slice(self.API.datacube, request.polytopes()) + time1 = time.time() + print("TIME TAKEN TO EXTRACT") + print(time1 - time0) + print(len(result.leaves)) + result.pprint() + + lats = [] + lons = [] + eccodes_lats = [] + eccodes_lons = [] + tol = 1e-8 + leaves = result.leaves + for i in range(len(leaves)): + cubepath = leaves[i].flatten() + lat = cubepath["latitude"][0] + lon = cubepath["longitude"][0] - 360 + lats.append(lat) + lons.append(lon) + + # # each variable in the netcdf file is a cube + # # cubes = iris.load('../../Downloads/votemper_ORAS5_1m_197902_grid_T_02.nc') + # # iris.save(cubes, '../../Downloads/votemper_ORAS5_1m_197902_grid_T_02.grib2') + # nearest_points = find_nearest_latlon("../../Downloads/votemper_ORAS5_1m_197902_grid_T_02.grib2", lat, lon) + # eccodes_lat = nearest_points[0][0]["lat"] + # eccodes_lon = nearest_points[0][0]["lon"] - 360 + # eccodes_lats.append(eccodes_lat) + # eccodes_lons.append(eccodes_lon) + # assert eccodes_lat - tol <= lat + # assert lat <= eccodes_lat + tol + # assert eccodes_lon - tol <= lon + # assert lon <= eccodes_lon + tol + + worldmap = gpd.read_file(gpd.datasets.get_path("naturalearth_lowres")) + fig, ax = plt.subplots(figsize=(12, 6)) + worldmap.plot(color="darkgrey", ax=ax) + + plt.scatter(lons, lats, s=18, c="red", cmap="YlOrRd") + plt.scatter(eccodes_lons, eccodes_lats, s=6, c="green") + plt.colorbar(label="Temperature") + plt.show() diff --git a/tests/test_override_md5_hash_options.py b/tests/test_override_md5_hash_options.py index 0233290af..237dc2ede 100644 --- a/tests/test_override_md5_hash_options.py +++ b/tests/test_override_md5_hash_options.py @@ -1,7 +1,6 @@ import pandas as pd import pytest -from polytope_feature.engine.hullslicer import HullSlicer from polytope_feature.polytope import Polytope, Request from polytope_feature.shapes import Box, Select, Span @@ -67,10 +66,8 @@ def test_fdb_datacube(self): ) self.fdbdatacube = gj.GribJump() - self.slicer = HullSlicer() self.API = Polytope( datacube=self.fdbdatacube, - engine=self.slicer, options=self.options, ) result = self.API.retrieve(request) diff --git a/tests/test_performance_icon_unstructured_fdb.py b/tests/test_performance_icon_unstructured_fdb.py new file mode 100644 index 000000000..0b79b52d7 --- /dev/null +++ b/tests/test_performance_icon_unstructured_fdb.py @@ -0,0 +1,93 @@ +import pygribjump as gj +import datetime +import math + +# import iris +import os +import time + +import geopandas as gpd +import matplotlib.pyplot as plt +import numpy as np +import pandas as pd +import pytest +import xarray as xr +from earthkit import data +from helper_functions import find_nearest_latlon + +from polytope_feature.polytope import Polytope, Request +from polytope_feature.shapes import Box, Point, Select, Polygon + +os.environ["FDB_HOME"] = "/Users/male/git/fdb-new-home" + + +engine_options = { + "step": "hullslicer", + "date": "hullslicer", + "levtype": "hullslicer", + "param": "hullslicer", + "latitude": "quadtree", + "longitude": "quadtree", +} +print("SETTING UP THE XARRAY") +time_now = time.time() + + +ds = data.from_source("file", "../../Downloads/icon_global_icosahedral_single-level_2025011000_000_T_2M.grib2") + +grid = xr.open_dataset("../../Downloads/icon_grid_0026_R03B07_G.nc", engine="netcdf4") + +print(time.time() - time_now) +arr = ds.to_xarray(engine="cfgrib").t2m + +longitudes = grid.clon.values * 180 / math.pi +latitudes = grid.clat.values * 180 / math.pi + +points = list(zip(latitudes, longitudes)) +print((min(latitudes), max(latitudes), min(longitudes), max(longitudes))) +print("FINISH SETTING UP POINTS") +options = { + "axis_config": [ + {"axis_name": "step", "transformations": [{"name": "type_change", "type": "int"}]}, + { + "axis_name": "date", + "transformations": [{"name": "merge", "other_axis": "time", "linkers": ["T", "00"]}], + }, + { + "axis_name": "values", + "transformations": [ + { + "name": "mapper", + "type": "irregular", + "resolution": 0, + "axes": ["latitude", "longitude"], + "md5_hash": "f68071a8ac9bae4e965822afb963c04f", + } + ], + }, + ], + "pre_path": {"date": "20250110"}, +} + + +tri_side = 80 +triangle = Polygon(["latitude", "longitude"], [[0, tri_side], [0, 0], [tri_side, 0]]) + +request = Request( + Select("date", [pd.Timestamp("20250110T0000")]), + Select("step", [0]), + Select("param", ["167"]), + Select("levtype", ["sfc"]), + Box(["latitude", "longitude"], [0, 0], [80, 80]), +) + +fdbdatacube = gj.GribJump() + +self_API = Polytope( + datacube=fdbdatacube, + options=options, + engine_options=engine_options, + point_cloud_options=points, +) +time0 = time.time() +result = self_API.retrieve(request) diff --git a/tests/test_point_nearest.py b/tests/test_point_nearest.py index b9fdc4bed..7b2e605b3 100644 --- a/tests/test_point_nearest.py +++ b/tests/test_point_nearest.py @@ -1,7 +1,6 @@ import pandas as pd import pytest -from polytope_feature.engine.hullslicer import HullSlicer from polytope_feature.polytope import Polytope, Request from polytope_feature.shapes import Point, Select @@ -60,10 +59,8 @@ def test_fdb_datacube(self): Point(["latitude", "longitude"], [[0.16, 0.176]], method="nearest"), ) self.fdbdatacube = gj.GribJump() - self.slicer = HullSlicer() self.API = Polytope( datacube=self.fdbdatacube, - engine=self.slicer, options=self.options, ) result = self.API.retrieve(request) @@ -86,10 +83,8 @@ def test_fdb_datacube_true_point(self): Point(["latitude", "longitude"], [[0.175746921078, 0.210608424337]], method="nearest"), ) self.fdbdatacube = gj.GribJump() - self.slicer = HullSlicer() self.API = Polytope( datacube=self.fdbdatacube, - engine=self.slicer, options=self.options, ) result = self.API.retrieve(request) @@ -113,10 +108,8 @@ def test_fdb_datacube_true_point_2(self): Point(["latitude", "longitude"], [[0.035149384216, 0.0]], method="nearest"), ) self.fdbdatacube = gj.GribJump() - self.slicer = HullSlicer() self.API = Polytope( datacube=self.fdbdatacube, - engine=self.slicer, options=self.options, ) result = self.API.retrieve(request) @@ -140,10 +133,8 @@ def test_fdb_datacube_true_point_3(self): Point(["latitude", "longitude"], [[0.035149384216, -0.01]], method="nearest"), ) self.fdbdatacube = gj.GribJump() - self.slicer = HullSlicer() self.API = Polytope( datacube=self.fdbdatacube, - engine=self.slicer, options=self.options, ) result = self.API.retrieve(request) @@ -169,10 +160,8 @@ def test_fdb_datacube_true_point_5(self): Point(["latitude", "longitude"], [[0.035149384216, 360 - 0.01]], method="nearest"), ) self.fdbdatacube = gj.GribJump() - self.slicer = HullSlicer() self.API = Polytope( datacube=self.fdbdatacube, - engine=self.slicer, options=self.options, ) result = self.API.retrieve(request) @@ -198,14 +187,11 @@ def test_fdb_datacube_true_point_4(self): Point(["latitude", "longitude"], [[0.035149384216, 359.97]], method="nearest"), ) self.fdbdatacube = gj.GribJump() - self.slicer = HullSlicer() self.API = Polytope( datacube=self.fdbdatacube, - engine=self.slicer, options=self.options, ) result = self.API.retrieve(request) - # result.pprint_2() assert len(result.leaves) == 1 assert result.leaves[0].values == (359.929906542056,) assert result.leaves[0].axis.name == "longitude" diff --git a/tests/test_point_shape.py b/tests/test_point_shape.py index d4d6eb98f..03a261217 100644 --- a/tests/test_point_shape.py +++ b/tests/test_point_shape.py @@ -2,7 +2,6 @@ import pandas as pd import xarray as xr -from polytope_feature.engine.hullslicer import HullSlicer from polytope_feature.polytope import Polytope, Request from polytope_feature.shapes import Point, Select @@ -19,9 +18,8 @@ def setup_method(self, method): "level": range(1, 130), }, ) - self.slicer = HullSlicer() options = {"compressed_axes_config": ["level", "step", "date"]} - self.API = Polytope(datacube=array, engine=self.slicer, options=options) + self.API = Polytope(datacube=array, options=options) def test_point(self): request = Request(Point(["step", "level"], [[3, 10]]), Select("date", ["2000-01-01"])) diff --git a/tests/test_polytope_extract.py b/tests/test_polytope_extract.py new file mode 100644 index 000000000..2bc7400fa --- /dev/null +++ b/tests/test_polytope_extract.py @@ -0,0 +1,51 @@ +import numpy as np +import xarray as xr + +from polytope_feature.polytope import Polytope, Request +from polytope_feature.shapes import Box + + +class TestPolytopeExtract: + def setup_method(self, method): + # Create a dataarray with 3 labelled axes using different index types + self.array = xr.DataArray( + np.random.randn(6, 129, 100), + dims=("step", "level", "values"), + coords={ + "step": [0, 3, 6, 9, 12, 15], + "level": range(1, 130), + "values": range(0, 100), + }, + ) + self.engine_options = { + "step": "hullslicer", + "level": "hullslicer", + "latitude": "quadtree", + "longitude": "quadtree", + } + self.quadtree_points = [[10, 10], [80, 10], [-5, 5], [5, 20], [5, 10], [50, 10]] + self.options = { + "axis_config": [ + { + "axis_name": "values", + "transformations": [ + {"name": "mapper", "type": "irregular", "resolution": 1280, "axes": ["latitude", "longitude"]} + ], + }, + ], + } + + # Testing different shapes + + def test_2D_box(self): + request = Request(Box(["step", "level"], [3, 10], [6, 11]), Box(["latitude", "longitude"], [0, 0], [20, 20])) + self.API = Polytope( + datacube=self.array, + options=self.options, + engine_options=self.engine_options, + point_cloud_options=self.quadtree_points, + ) + result = self.API.retrieve(request) + + result.pprint() + assert len(result.leaves) == 12 diff --git a/tests/test_polytope_extract_fdb.py b/tests/test_polytope_extract_fdb.py new file mode 100644 index 000000000..d466581b6 --- /dev/null +++ b/tests/test_polytope_extract_fdb.py @@ -0,0 +1,95 @@ +import pandas as pd +import pytest + +from polytope_feature.polytope import Polytope, Request +from polytope_feature.shapes import Box, Select + + +class TestPolytopeExtract: + def setup_method(self, method): + # from polytope.datacube.backends.fdb import FDBDatacube + + # Create a dataarray with 3 labelled axes using different index types + self.engine_options = { + "step": "hullslicer", + "levtype": "hullslicer", + "latitude": "quadtree", + "longitude": "quadtree", + "class": "hullslicer", + "date": "hullslicer", + "type": "hullslicer", + "stream": "hullslicer", + "param": "hullslicer", + "expver": "hullslicer", + "domain": "hullslicer", + } + self.quadtree_points = [[10, 10], [0.035149384216, 0.0], [80, 10], [-5, 5], [5, 20], [5, 10], [50, 10]] + self.options = { + "axis_config": [ + {"axis_name": "step", "transformations": [{"name": "type_change", "type": "int"}]}, + {"axis_name": "number", "transformations": [{"name": "type_change", "type": "int"}]}, + { + "axis_name": "date", + "transformations": [{"name": "merge", "other_axis": "time", "linkers": ["T", "00"]}], + }, + { + "axis_name": "values", + "transformations": [ + { + "name": "mapper", + "type": "irregular", + "axes": ["latitude", "longitude"], + "md5_hash": "158db321ae8e773681eeb40e0a3d350f", + } + ], + }, + ], + "compressed_axes_config": [ + "longitude", + "latitude", + "levtype", + "step", + "date", + "domain", + "expver", + "param", + "class", + "stream", + "type", + ], + "pre_path": {"class": "od", "expver": "0001", "levtype": "sfc", "stream": "oper"}, + } + + # Testing different shapes + @pytest.mark.fdb + def test_2D_box(self): + import pygribjump as gj + + request = Request( + Select("step", [0]), + Select("levtype", ["sfc"]), + Select("date", [pd.Timestamp("20230625T120000")]), + Select("domain", ["g"]), + Select("expver", ["0001"]), + Select("param", ["167"]), + Select("class", ["od"]), + Select("stream", ["oper"]), + Select("type", ["an"]), + Box(["latitude", "longitude"], [0, -0.1], [10, 10]), + ) + self.fdbdatacube = gj.GribJump() + self.API = Polytope( + datacube=self.fdbdatacube, + options=self.options, + engine_options=self.engine_options, + point_cloud_options=self.quadtree_points, + ) + result = self.API.retrieve(request) + + assert len(result.leaves) == 3 + assert result.leaves[0].flatten()["longitude"] == (0,) + assert result.leaves[0].flatten()["latitude"] == (0.035149384216,) + assert result.leaves[1].flatten()["longitude"] == (10,) + assert result.leaves[1].flatten()["latitude"] == (5,) + assert result.leaves[2].flatten()["longitude"] == (10,) + assert result.leaves[2].flatten()["latitude"] == (10,) diff --git a/tests/test_quad_tree.py b/tests/test_quad_tree.py new file mode 100644 index 000000000..d84d52d5f --- /dev/null +++ b/tests/test_quad_tree.py @@ -0,0 +1,162 @@ +import pytest + +from polytope_feature.datacube.quad_tree import QuadNode +from polytope_feature.datacube.quadtree_additional_operations import query_polygon +from polytope_feature.engine.quadtree_slicer import QuadTreeSlicer +from polytope_feature.engine.slicing_tools import slice_in_two +from polytope_feature.polytope import Polytope +from polytope_feature.shapes import Box, ConvexPolytope + + +class TestQuadTreeSlicer: + def setup_method(self, method): + import pygribjump as gj + + self.options = { + "axis_config": [ + {"axis_name": "step", "transformations": [{"name": "type_change", "type": "int"}]}, + {"axis_name": "number", "transformations": [{"name": "type_change", "type": "int"}]}, + { + "axis_name": "date", + "transformations": [{"name": "merge", "other_axis": "time", "linkers": ["T", "00"]}], + }, + { + "axis_name": "values", + "transformations": [ + {"name": "mapper", "type": "irregular", "resolution": 1280, "axes": ["latitude", "longitude"]} + ], + }, + ], + "compressed_axes_config": [ + "longitude", + "latitude", + "levtype", + "step", + "date", + "domain", + "expver", + "param", + "class", + "stream", + "type", + ], + "pre_path": {"class": "od", "expver": "0001", "levtype": "sfc", "stream": "oper"}, + } + self.fdbdatacube = gj.GribJump() + + @pytest.mark.fdb + def test_quad_tree_slicer(self): + points = [[10, 10], [80, 10], [-5, 5], [5, 20], [5, 10], [50, 10]] + slicer = QuadTreeSlicer(points) + # slicer.quad_tree.pprint() + pass + + @pytest.mark.fdb + def test_quad_tree_query_polygon(self): + points = [[10, 10], [80, 10], [-5, 5], [5, 20], [5, 10], [50, 10]] + slicer = QuadTreeSlicer(points) + polytope = Box(["lat", "lon"], [1, 1], [20, 30]).polytope()[0] + # results = slicer.quad_tree.query_polygon(polytope) + results = query_polygon(points, slicer.quad_tree, 0, polytope, results=None) + assert len(results) == 3 + assert (10, 10) in [slicer.points[node] for node in results] + assert (5, 10) in [slicer.points[node] for node in results] + assert (5, 20) in [slicer.points[node] for node in results] + points = [[10, 10], [80, 10], [-5, 5], [5, 50], [5, 10], [50, 10], [2, 10], [15, 15]] + slicer = QuadTreeSlicer(points) + polytope = ConvexPolytope(["lat", "lon"], [[-10, 1], [20, 1], [5, 20]]) + # results = slicer.quad_tree.query_polygon(polytope) + results = query_polygon(points, slicer.quad_tree, 0, polytope, results=None) + assert len(results) == 4 + assert (-5, 5) in [slicer.points[node] for node in results] + assert (5, 10) in [slicer.points[node] for node in results] + assert (10, 10) in [slicer.points[node] for node in results] + assert (2, 10) in [slicer.points[node] for node in results] + + @pytest.mark.fdb + def test_slice_in_two_vertically(self): + polytope = Box(["lat", "lon"], [0, 0], [2, 2]).polytope()[0] + lower, upper = slice_in_two(polytope, 1, 0) + assert lower.points == [[0, 0], [1.0, 0.0], [1.0, 2.0], [0, 2]] + assert upper.points == [[1.0, 0.0], [2, 0], [2, 2], [1.0, 2.0]] + + @pytest.mark.fdb + def test_slice_in_two_horizontally(self): + polytope = Box(["lat", "lon"], [0, 0], [2, 2]).polytope()[0] + lower, upper = slice_in_two(polytope, 1, 1) + assert lower.points == [[0, 0], [2, 0], [2.0, 1.0], [0.0, 1.0]] + assert upper.points == [[2, 2], [0, 2], [0.0, 1.0], [2.0, 1.0]] + + @pytest.mark.fdb + def test_quad_node_is_contained_in_box(self): + node = QuadNode([1, 1], 0) + polytope = Box(["lat", "lon"], [0, 0], [2, 2]).polytope()[0] + assert node.is_contained_in(polytope) + second_node = QuadNode([3, 3], 0) + assert not second_node.is_contained_in(polytope) + third_node = QuadNode([1, 0], 0) + assert third_node.is_contained_in(polytope) + + @pytest.mark.fdb + def test_quad_node_is_contained_in_triangle(self): + node = QuadNode([1, 1], 0) + polytope = ConvexPolytope(["lat", "lon"], [[0, 0], [1, 1], [2, 0]]) + assert node.is_contained_in(polytope) + node = QuadNode([1, 0.5], 0) + assert node.is_contained_in(polytope) + second_node = QuadNode([3, 3], 0) + assert not second_node.is_contained_in(polytope) + third_node = QuadNode([1, 0], 0) + assert third_node.is_contained_in(polytope) + third_node = QuadNode([0.1, 0.5], 0) + assert not third_node.is_contained_in(polytope) + + @pytest.mark.fdb + def test_quad_tree_slicer_extract(self): + points = [[10, 10], [80, 10], [-5, 5], [5, 20], [5, 10], [50, 10]] + polytope = Box(["latitude", "longitude"], [1, 1], [20, 30]).polytope()[0] + self.API = Polytope( + datacube=self.fdbdatacube, + options=self.options, + engine_options={"latitude": "quadtree", "longitude": "quadtree"}, + point_cloud_options=points, + ) + tree = self.API.engines["quadtree"].extract(self.API.datacube, [polytope]) + assert len(tree.leaves) == 3 + # tree.pprint() + points = [[10, 10], [80, 10], [-5, 5], [5, 50], [5, 10], [50, 10], [2, 10], [15, 15]] + polytope = ConvexPolytope(["latitude", "longitude"], [[-10, 1], [20, 1], [5, 20]]) + tree = self.API.engines["quadtree"].extract(self.API.datacube, [polytope]) + assert len(tree.leaves) == 4 + # tree.pprint() + + # @pytest.mark.skip("performance test") + @pytest.mark.fdb + def test_large_scale_extraction(self): + import time + + import numpy as np + + x = np.linspace(0, 100, 1000) + y = np.linspace(0, 100, 1000) + # create the mesh based on these arrays + X, Y = np.meshgrid(x, y) + X = X.reshape((np.prod(X.shape),)) + Y = Y.reshape((np.prod(Y.shape),)) + coords = zip(X, Y) + points = [list(coord) for coord in coords] + time0 = time.time() + polytope = Box(["latitude", "longitude"], [1, 1], [20, 30]).polytope()[0] + self.API = Polytope( + datacube=self.fdbdatacube, + options=self.options, + engine_options={"latitude": "quadtree", "longitude": "quadtree"}, + point_cloud_options=points, + ) + print(time.time() - time0) + time1 = time.time() + tree = self.API.engines["quadtree"].extract(self.API.datacube, [polytope]) + print(time.time() - time1) # = 5.919436931610107 + print(len(tree.leaves)) # = 55100 + # NOTE: maybe for 2D qhull here, scipy is not the fastest + # but use shewchuk's triangle algo: https://www.cs.cmu.edu/~quake/triangle.html? diff --git a/tests/test_quadtree_edge_cases.py b/tests/test_quadtree_edge_cases.py new file mode 100644 index 000000000..237ceac4c --- /dev/null +++ b/tests/test_quadtree_edge_cases.py @@ -0,0 +1,56 @@ +import pytest + +from polytope_feature.polytope import Polytope +from polytope_feature.shapes import Box + + +class TestQuadTreeSlicer: + def setup_method(self, method): + import pygribjump as gj + + self.options = { + "axis_config": [ + {"axis_name": "step", "transformations": [{"name": "type_change", "type": "int"}]}, + {"axis_name": "number", "transformations": [{"name": "type_change", "type": "int"}]}, + { + "axis_name": "date", + "transformations": [{"name": "merge", "other_axis": "time", "linkers": ["T", "00"]}], + }, + { + "axis_name": "values", + "transformations": [ + {"name": "mapper", "type": "irregular", "resolution": 1280, "axes": ["latitude", "longitude"]} + ], + }, + ], + "compressed_axes_config": [ + "longitude", + "latitude", + "levtype", + "step", + "date", + "domain", + "expver", + "param", + "class", + "stream", + "type", + ], + "pre_path": {"class": "od", "expver": "0001", "levtype": "sfc", "stream": "oper"}, + } + self.fdbdatacube = gj.GribJump() + + @pytest.mark.fdb + def test_quad_tree_slicer_extract(self): + points = [[10, 10], [80, 10], [-5, 5], [5, 20], [5, 10], [50, 10], [0.035149384216, 0.0]] + polytope = Box(["latitude", "longitude"], [0, 0], [15, 15]).polytope()[0] + self.API = Polytope( + datacube=self.fdbdatacube, + options=self.options, + engine_options={"latitude": "quadtree", "longitude": "quadtree"}, + point_cloud_options=points, + ) + tree = self.API.engines["quadtree"].extract(self.API.datacube, [polytope]) + tree.pprint() + assert len(tree.leaves) == 3 + assert set([leaf.flatten()["values"] for leaf in tree.leaves]) == set([(0,), (4,), (6,)]) diff --git a/tests/test_quadtree_indices.py b/tests/test_quadtree_indices.py new file mode 100644 index 000000000..c84c4ee4d --- /dev/null +++ b/tests/test_quadtree_indices.py @@ -0,0 +1,68 @@ +import pytest + +from polytope_feature.polytope import Polytope +from polytope_feature.shapes import Box, ConvexPolytope + + +class TestQuadTreeSlicer: + def setup_method(self, method): + import pygribjump as gj + + self.options = { + "axis_config": [ + {"axis_name": "step", "transformations": [{"name": "type_change", "type": "int"}]}, + {"axis_name": "number", "transformations": [{"name": "type_change", "type": "int"}]}, + { + "axis_name": "date", + "transformations": [{"name": "merge", "other_axis": "time", "linkers": ["T", "00"]}], + }, + { + "axis_name": "values", + "transformations": [ + {"name": "mapper", "type": "irregular", "resolution": 1280, "axes": ["latitude", "longitude"]} + ], + }, + ], + "compressed_axes_config": [ + "longitude", + "latitude", + "levtype", + "step", + "date", + "domain", + "expver", + "param", + "class", + "stream", + "type", + ], + "pre_path": {"class": "od", "expver": "0001", "levtype": "sfc", "stream": "oper"}, + } + self.fdbdatacube = gj.GribJump() + + @pytest.mark.fdb + def test_quad_tree_slicer_extract(self): + points = [[10, 10], [80, 10], [-5, 5], [5, 20], [5, 10], [50, 10]] + polytope = Box(["latitude", "longitude"], [1, 1], [20, 30]).polytope()[0] + self.API = Polytope( + datacube=self.fdbdatacube, + options=self.options, + engine_options={"latitude": "quadtree", "longitude": "quadtree"}, + point_cloud_options=points, + ) + tree = self.API.engines["quadtree"].extract(self.API.datacube, [polytope]) + assert len(tree.leaves) == 3 + assert set([leaf.flatten()["values"] for leaf in tree.leaves]) == set([(0,), (3,), (4,)]) + tree.pprint() + points = [[10, 10], [80, 10], [-5, 5], [5, 50], [5, 10], [50, 10], [2, 10], [15, 15]] + polytope = ConvexPolytope(["latitude", "longitude"], [[-10, 1], [20, 1], [5, 20]]) + self.API = Polytope( + datacube=self.fdbdatacube, + options=self.options, + engine_options={"latitude": "quadtree", "longitude": "quadtree"}, + point_cloud_options=points, + ) + tree = self.API.engines["quadtree"].extract(self.API.datacube, [polytope]) + assert len(tree.leaves) == 4 + assert set([leaf.flatten()["values"] for leaf in tree.leaves]) == set([(0,), (2,), (4,), (6,)]) + tree.pprint() diff --git a/tests/test_quadtree_optimisation.py b/tests/test_quadtree_optimisation.py new file mode 100644 index 000000000..81a57e926 --- /dev/null +++ b/tests/test_quadtree_optimisation.py @@ -0,0 +1,63 @@ +import pytest + +from polytope_feature.datacube.quadtree_additional_operations import query_polygon +from polytope_feature.engine.quadtree_slicer import QuadTreeSlicer +from polytope_feature.shapes import Box + + +class TestQuadTreeSlicer: + def setup_method(self, method): + import pygribjump as gj + + self.options = { + "axis_config": [ + {"axis_name": "step", "transformations": [{"name": "type_change", "type": "int"}]}, + {"axis_name": "number", "transformations": [{"name": "type_change", "type": "int"}]}, + { + "axis_name": "date", + "transformations": [{"name": "merge", "other_axis": "time", "linkers": ["T", "00"]}], + }, + { + "axis_name": "values", + "transformations": [ + {"name": "mapper", "type": "irregular", "resolution": 1280, "axes": ["latitude", "longitude"]} + ], + }, + ], + "compressed_axes_config": [ + "longitude", + "latitude", + "levtype", + "step", + "date", + "domain", + "expver", + "param", + "class", + "stream", + "type", + ], + "pre_path": {"class": "od", "expver": "0001", "levtype": "sfc", "stream": "oper"}, + } + self.fdbdatacube = gj.GribJump() + + @pytest.mark.fdb + def test_quad_tree_slicer(self): + points = [[10, 10], [80, 10], [-5, 5], [5, 20], [5, 10], [50, 10]] + slicer = QuadTreeSlicer(points) + # slicer.quad_tree.pprint() + pass + + @pytest.mark.fdb + def test_quad_tree_query_polygon(self): + points = [[10, 10], [80, 10], [-5, 5], [5, 20], [5, 10], [50, 10]] + slicer = QuadTreeSlicer(points) + polytope = Box(["lat", "lon"], [0, 0], [90, 45]).polytope()[0] + # results = slicer.quad_tree.query_polygon(polytope) + results = query_polygon(points, slicer.quad_tree, 0, polytope, results=None) + assert len(results) == 5 + assert (10, 10) in [slicer.points[node] for node in results] + assert (5, 10) in [slicer.points[node] for node in results] + assert (5, 20) in [slicer.points[node] for node in results] + assert (80, 10) in [slicer.points[node] for node in results] + assert (50, 10) in [slicer.points[node] for node in results] diff --git a/tests/test_qubed_extraction.py b/tests/test_qubed_extraction.py new file mode 100644 index 000000000..7e72975d3 --- /dev/null +++ b/tests/test_qubed_extraction.py @@ -0,0 +1,67 @@ +from qubed import Qube +import requests +from polytope_feature.datacube.datacube_axis import PandasTimedeltaDatacubeAxis, PandasTimestampDatacubeAxis, UnsliceableDatacubeAxis + +from polytope_feature.shapes import ConvexPolytope + + +fdb_tree = Qube.from_json(requests.get( + "https://github.com/ecmwf/qubed/raw/refs/heads/main/tests/example_qubes/climate_dt.json").json()) + +print(fdb_tree.axes().keys()) + + +combi_polytopes = [ + ConvexPolytope(["param"], [["168"]]), + ConvexPolytope(["time"], [["0000"], ["1200"]]), + ConvexPolytope(["resolution"], [["high"]]), + ConvexPolytope(["type"], [["fc"]]), + ConvexPolytope(["model"], ['ifs-nemo']), + ConvexPolytope(["stream"], [["clte"]]), + ConvexPolytope(["realization"], ["1"]), + ConvexPolytope(["expver"], [['0001']]), + ConvexPolytope(["experiment"], [['ssp3-7.0']]), + ConvexPolytope(["generation"], [["1"]]), + ConvexPolytope(["levtype"], [["sfc"]]), + ConvexPolytope(["activity"], [["scenariomip"]]), + ConvexPolytope(["dataset"], [["climate-dt"]]), + ConvexPolytope(["class"], [["d1"]]), + ConvexPolytope(["date"], [["20190221", "20190223"]]) +] + +datacube_axes = {"param": UnsliceableDatacubeAxis(), + "time": PandasTimedeltaDatacubeAxis(), + "resolution": UnsliceableDatacubeAxis(), + "type": UnsliceableDatacubeAxis(), + "model": UnsliceableDatacubeAxis(), + "stream": UnsliceableDatacubeAxis(), + "realization": UnsliceableDatacubeAxis(), + "expver": UnsliceableDatacubeAxis(), + "experiment": UnsliceableDatacubeAxis(), + "generation": UnsliceableDatacubeAxis(), + "levtype": UnsliceableDatacubeAxis(), + "activity": UnsliceableDatacubeAxis(), + "dataset": UnsliceableDatacubeAxis(), + "class": UnsliceableDatacubeAxis(), + "date": PandasTimestampDatacubeAxis()} + + +# TODO: treat the transformations to talk to the qubed tree, maybe do it + + +# TODO: start iterating fdb_tree and creating a new request tree + + +# print(fdb_tree.) + + +# Select("step", [0]), +# Select("levtype", ["sfc"]), +# Select("date", [pd.Timestamp("20231102T000000")]), +# Select("domain", ["g"]), +# Select("expver", ["0001"]), +# Select("param", ["167"]), +# Select("class", ["od"]), +# Select("stream", ["oper"]), +# Select("type", ["fc"]), +# Box(["latitude", "longitude"], [0, 0], [80, 80]), diff --git a/tests/test_regular_grid.py b/tests/test_regular_grid.py index 7eaaeccb5..bea22f265 100644 --- a/tests/test_regular_grid.py +++ b/tests/test_regular_grid.py @@ -2,7 +2,6 @@ import pytest from helper_functions import download_test_data, find_nearest_latlon -from polytope_feature.engine.hullslicer import HullSlicer from polytope_feature.polytope import Polytope, Request from polytope_feature.shapes import Disk, Select @@ -76,10 +75,8 @@ def test_regular_grid(self): Select("number", ["0", "1"]), ) self.fdbdatacube = gj.GribJump() - self.slicer = HullSlicer() self.API = Polytope( datacube=self.fdbdatacube, - engine=self.slicer, options=self.options, ) result = self.API.retrieve(request) diff --git a/tests/test_request_tree.py b/tests/test_request_tree.py index 038599724..d0976b941 100644 --- a/tests/test_request_tree.py +++ b/tests/test_request_tree.py @@ -274,7 +274,8 @@ def test_get_ancestors(self): child1.add_child(grandchild1) root_node1 = TensorIndexTree() root_node1.add_child(child1) - assert greatgrandchild1.get_ancestors() == SortedList([greatgrandchild1, grandchild1, child1]) + # assert greatgrandchild1.get_ancestors() == SortedList([greatgrandchild1, grandchild1, child1]) + assert greatgrandchild1.get_ancestors() == SortedList([child1, grandchild1, greatgrandchild1]) def test_add_or_get_child(self): axis1 = IntDatacubeAxis() @@ -290,7 +291,7 @@ def test_add_or_get_child(self): root_node = TensorIndexTree() root_node.add_child(child1) assert root_node.create_child(axis1, 0, [])[0] == child1 - assert root_node.create_child(axis2, (), [])[0].parent == root_node + assert root_node.create_child(axis2, 0, [])[0].parent == root_node def test_eq(self): axis1 = IntDatacubeAxis() diff --git a/tests/test_request_trees_after_slicing.py b/tests/test_request_trees_after_slicing.py index 2aa28b4a7..af13ab352 100644 --- a/tests/test_request_trees_after_slicing.py +++ b/tests/test_request_trees_after_slicing.py @@ -3,7 +3,6 @@ from polytope_feature.datacube.backends.xarray import XArrayDatacube from polytope_feature.datacube.datacube_axis import IntDatacubeAxis -from polytope_feature.engine.hullslicer import HullSlicer from polytope_feature.polytope import Polytope from polytope_feature.shapes import Box @@ -19,14 +18,13 @@ def setup_method(self, method): }, ) self.xarraydatacube = XArrayDatacube(array) - self.slicer = HullSlicer() options = {"compressed_axes_config": ["level", "step"]} - self.API = Polytope(datacube=array, engine=self.slicer, options=options) + self.API = Polytope(datacube=array, options=options) def test_path_values(self): box = Box(["step", "level"], [3.0, 1.0], [6.0, 3.0]) polytope = box.polytope() - request = self.slicer.extract(self.xarraydatacube, polytope) + request = self.API.slice(self.xarraydatacube, polytope) datacube_path = request.leaves[0].flatten() request.pprint() assert datacube_path.values() == tuple([tuple([3.0]), tuple([1.0, 2, 3])]) @@ -35,7 +33,7 @@ def test_path_values(self): def test_path_keys(self): box = Box(["step", "level"], [3.0, 1.0], [6.0, 3.0]) polytope = box.polytope() - request = self.slicer.extract(self.xarraydatacube, polytope) + request = self.API.slice(self.xarraydatacube, polytope) datacube_path = request.leaves[0].flatten() assert datacube_path.keys()[0] == "step" assert datacube_path.keys()[1] == "level" @@ -43,14 +41,14 @@ def test_path_keys(self): def test_path_pprint(self): box = Box(["step", "level"], [3.0, 1.0], [6.0, 3.0]) polytope = box.polytope() - request = self.slicer.extract(self.xarraydatacube, polytope) + request = self.API.slice(self.xarraydatacube, polytope) datacube_path = request.leaves[0].flatten() datacube_path.pprint() def test_flatten(self): box = Box(["step", "level"], [3.0, 1.0], [6.0, 3.0]) polytope = box.polytope() - request = self.slicer.extract(self.xarraydatacube, polytope) + request = self.API.slice(self.xarraydatacube, polytope) path = request.leaves[0].flatten() request.pprint() assert path["step"] == tuple([3.0]) @@ -59,7 +57,7 @@ def test_flatten(self): def test_add_child(self): box = Box(["step", "level"], [3.0, 1.0], [6.0, 3.0]) polytope = box.polytope() - request = self.slicer.extract(self.xarraydatacube, polytope) + request = self.API.slice(self.xarraydatacube, polytope) request1 = request.leaves[0] request2 = request.leaves[0] # Test adding child @@ -77,13 +75,13 @@ def test_add_child(self): def test_pprint(self): box = Box(["step", "level"], [3.0, 1.0], [6.0, 3.0]) polytope = box.polytope() - request = self.slicer.extract(self.xarraydatacube, polytope) + request = self.API.slice(self.xarraydatacube, polytope) request.pprint() def test_remove_branch(self): box = Box(["step", "level"], [3.0, 1.0], [6.0, 3.0]) polytope = box.polytope() - request = self.slicer.extract(self.xarraydatacube, polytope) + request = self.API.slice(self.xarraydatacube, polytope) prev_request_size = len(request.leaves) request.leaves[0].remove_branch() new_request_size = len(request.leaves) diff --git a/tests/test_reverse_transformation.py b/tests/test_reverse_transformation.py index f1d9b24ed..1b752ae3a 100644 --- a/tests/test_reverse_transformation.py +++ b/tests/test_reverse_transformation.py @@ -21,7 +21,7 @@ def setup_method(self, method): "compressed_axes_config": ["lat"], } self.slicer = HullSlicer() - self.API = Polytope(datacube=array, engine=self.slicer, options=options) + self.API = Polytope(datacube=array, options=options) def test_reverse_transformation(self): request = Request(Select("lat", [1, 2, 3])) diff --git a/tests/test_shapes.py b/tests/test_shapes.py index ac730905a..b6408ccc8 100644 --- a/tests/test_shapes.py +++ b/tests/test_shapes.py @@ -3,7 +3,6 @@ import pytest import xarray as xr -from polytope_feature.engine.hullslicer import HullSlicer from polytope_feature.polytope import Polytope, Request from polytope_feature.shapes import All, Select, Span @@ -25,10 +24,8 @@ def setup_method(self, method): "axis_config": [{"axis_name": "longitude", "transformations": [{"name": "cyclic", "range": [0, 360]}]}], "compressed_axes_config": ["date", "step", "level", "longitude"], } - self.slicer = HullSlicer() self.API = Polytope( datacube=array, - engine=self.slicer, options=self.options, ) @@ -76,7 +73,6 @@ def test_all_mapper_cyclic(self): "pre_path": {"class": "od", "expver": "0001", "levtype": "sfc", "stream": "oper"}, } self.fdbdatacube = gj.GribJump() - self.slicer = HullSlicer() request = Request( Select("step", [11]), @@ -91,7 +87,7 @@ def test_all_mapper_cyclic(self): Span("latitude", 89.9, 90), All("longitude"), ) - self.API = Polytope(datacube=self.fdbdatacube, engine=self.slicer, options=self.options) + self.API = Polytope(datacube=self.fdbdatacube, options=self.options) result = self.API.retrieve(request) result.pprint() assert len(result.leaves) == 1 diff --git a/tests/test_slice_date_range_fdb.py b/tests/test_slice_date_range_fdb.py index dae81cb22..ed90b903c 100644 --- a/tests/test_slice_date_range_fdb.py +++ b/tests/test_slice_date_range_fdb.py @@ -1,7 +1,6 @@ import pandas as pd import pytest -from polytope_feature.engine.hullslicer import HullSlicer from polytope_feature.polytope import Polytope, Request from polytope_feature.shapes import Box, Disk, Select, Span @@ -61,10 +60,8 @@ def test_fdb_datacube(self): ) self.fdbdatacube = gj.GribJump() - self.slicer = HullSlicer() self.API = Polytope( datacube=self.fdbdatacube, - engine=self.slicer, options=self.options, ) result = self.API.retrieve(request) @@ -91,10 +88,8 @@ def test_fdb_datacube_select_non_existing_last(self): ) self.fdbdatacube = gj.GribJump() - self.slicer = HullSlicer() self.API = Polytope( datacube=self.fdbdatacube, - engine=self.slicer, options=self.options, ) result = self.API.retrieve(request) @@ -121,10 +116,8 @@ def test_fdb_datacube_select_non_existing_first(self): ) self.fdbdatacube = gj.GribJump() - self.slicer = HullSlicer() self.API = Polytope( datacube=self.fdbdatacube, - engine=self.slicer, options=self.options, ) result = self.API.retrieve(request) @@ -134,37 +127,37 @@ def test_fdb_datacube_select_non_existing_first(self): assert len(result.leaves[i].result) == 3 @pytest.mark.fdb - def test_fdb_datacube_select_completely_non_existing(self): + def test_fdb_datacube_disk(self): import pygribjump as gj request = Request( Select("step", [0]), Select("levtype", ["sfc"]), - Select("date", [pd.Timestamp("20230624T120000"), pd.Timestamp("20230626T120000")]), + Span("date", pd.Timestamp("20230625T120000"), pd.Timestamp("20230626T120000")), Select("domain", ["g"]), Select("expver", ["0001"]), Select("param", ["167"]), Select("class", ["od"]), Select("stream", ["oper"]), Select("type", ["an"]), - Box(["latitude", "longitude"], [0, 0], [0.2, 0.2]), + Disk(["latitude", "longitude"], [0, 0], [0.1, 0.1]), ) self.fdbdatacube = gj.GribJump() - self.slicer = HullSlicer() self.API = Polytope( datacube=self.fdbdatacube, - engine=self.slicer, options=self.options, ) result = self.API.retrieve(request) result.pprint() - assert len(result.leaves) == 1 - for i in range(len(result.leaves)): - assert len(result.leaves[i].result) == 0 + assert len(result.leaves) == 2 + assert len(result.leaves[0].result) == 3 + assert len(result.leaves[1].result) == 3 + assert len(result.leaves[0].values) == 3 + assert len(result.leaves[1].values) == 3 @pytest.mark.fdb - def test_fdb_datacube_disk(self): + def test_fdb_datacube_disk_2(self): import pygribjump as gj request = Request( @@ -177,54 +170,48 @@ def test_fdb_datacube_disk(self): Select("class", ["od"]), Select("stream", ["oper"]), Select("type", ["an"]), - Disk(["latitude", "longitude"], [0, 0], [0.1, 0.1]), + Disk(["latitude", "longitude"], [0.05, 0.070148090413], [0.1, 0.15]), ) self.fdbdatacube = gj.GribJump() - self.slicer = HullSlicer() self.API = Polytope( datacube=self.fdbdatacube, - engine=self.slicer, options=self.options, ) result = self.API.retrieve(request) result.pprint() - assert len(result.leaves) == 2 + assert len(result.leaves) == 3 assert len(result.leaves[0].result) == 3 - assert len(result.leaves[1].result) == 3 + assert len(result.leaves[1].result) == 5 + assert len(result.leaves[2].result) == 3 assert len(result.leaves[0].values) == 3 - assert len(result.leaves[1].values) == 3 + assert len(result.leaves[1].values) == 5 + assert len(result.leaves[2].values) == 3 @pytest.mark.fdb - def test_fdb_datacube_disk_2(self): + def test_fdb_datacube_select_completely_non_existing(self): import pygribjump as gj request = Request( Select("step", [0]), Select("levtype", ["sfc"]), - Span("date", pd.Timestamp("20230625T120000"), pd.Timestamp("20230626T120000")), + Select("date", [pd.Timestamp("20230624T120000"), pd.Timestamp("20230626T120000")]), Select("domain", ["g"]), Select("expver", ["0001"]), Select("param", ["167"]), Select("class", ["od"]), Select("stream", ["oper"]), Select("type", ["an"]), - Disk(["latitude", "longitude"], [0.05, 0.070148090413], [0.1, 0.15]), + Box(["latitude", "longitude"], [0, 0], [0.2, 0.2]), ) self.fdbdatacube = gj.GribJump() - self.slicer = HullSlicer() self.API = Polytope( datacube=self.fdbdatacube, - engine=self.slicer, options=self.options, ) result = self.API.retrieve(request) result.pprint() - assert len(result.leaves) == 3 - assert len(result.leaves[0].result) == 3 - assert len(result.leaves[1].result) == 5 - assert len(result.leaves[2].result) == 3 - assert len(result.leaves[0].values) == 3 - assert len(result.leaves[1].values) == 5 - assert len(result.leaves[2].values) == 3 + assert len(result.leaves) == 1 + for i in range(len(result.leaves)): + assert len(result.leaves[i].result) == 0 diff --git a/tests/test_slice_date_range_fdb_v2.py b/tests/test_slice_date_range_fdb_v2.py index fcacf6b7c..cb76be722 100644 --- a/tests/test_slice_date_range_fdb_v2.py +++ b/tests/test_slice_date_range_fdb_v2.py @@ -2,7 +2,6 @@ import pandas as pd import pytest -from polytope_feature.engine.hullslicer import HullSlicer from polytope_feature.polytope import Polytope, Request from polytope_feature.shapes import Select, Span @@ -67,10 +66,8 @@ def test_fdb_datacube(self): Select("number", ["0"]), ) self.fdbdatacube = gj.GribJump() - self.slicer = HullSlicer() self.API = Polytope( datacube=self.fdbdatacube, - engine=self.slicer, options=self.options, ) result = self.API.retrieve(request) diff --git a/tests/test_slicer_engine.py b/tests/test_slicer_engine.py index 289eeab6b..ba04b2672 100644 --- a/tests/test_slicer_engine.py +++ b/tests/test_slicer_engine.py @@ -3,7 +3,6 @@ from polytope_feature.datacube.backends.xarray import XArrayDatacube from polytope_feature.datacube.tensor_index_tree import TensorIndexTree -from polytope_feature.engine.hullslicer import HullSlicer from polytope_feature.polytope import Polytope from polytope_feature.shapes import Box @@ -19,14 +18,13 @@ def setup_method(self, method): }, ) self.xarraydatacube = XArrayDatacube(array) - self.slicer = HullSlicer() options = {"compressed_axes_config": ["level", "step"]} - self.API = Polytope(datacube=array, engine=self.slicer, options=options) + self.API = Polytope(datacube=array, options=options) def test_extract(self): box = Box(["step", "level"], [3.0, 1.0], [6.0, 3.0]) polytope = box.polytope() - request = self.slicer.extract(self.xarraydatacube, polytope) + request = self.API.slice(self.xarraydatacube, polytope) assert request.axis == TensorIndexTree.root assert request.parent is None assert request.values is tuple() diff --git a/tests/test_slicer_era5.py b/tests/test_slicer_era5.py index 9f11bcef8..c3e450c07 100644 --- a/tests/test_slicer_era5.py +++ b/tests/test_slicer_era5.py @@ -3,7 +3,6 @@ from earthkit import data from helper_functions import download_test_data -from polytope_feature.engine.hullslicer import HullSlicer from polytope_feature.polytope import Polytope, Request from polytope_feature.shapes import Box, Select @@ -15,14 +14,12 @@ def setup_method(self, method): ds = data.from_source("file", "./tests/data/era5-levels-members.grib") array = ds.to_xarray(engine="cfgrib").isel(step=0).t - self.slicer = HullSlicer() options = { "axis_config": [{"axis_name": "latitude", "transformations": [{"name": "reverse", "is_reverse": True}]}], "compressed_axes_config": ["number", "time", "latitude", "longitude", "step", "isobaricInhPa"], } self.API = Polytope( datacube=array, - engine=self.slicer, options=options, ) diff --git a/tests/test_slicer_xarray.py b/tests/test_slicer_xarray.py index b36cd051c..3aa8d3f7c 100644 --- a/tests/test_slicer_xarray.py +++ b/tests/test_slicer_xarray.py @@ -2,7 +2,6 @@ import pandas as pd import xarray as xr -from polytope_feature.engine.hullslicer import HullSlicer from polytope_feature.polytope import Polytope, Request from polytope_feature.shapes import Box, Select, Span @@ -19,9 +18,8 @@ def setup_method(self, method): "level": range(1, 130), }, ) - self.slicer = HullSlicer() options = {"compressed_axes_config": ["date", "step", "level"]} - self.API = Polytope(datacube=array, engine=self.slicer, options=options) + self.API = Polytope(datacube=array, options=options) def test_2D_box(self): request = Request(Box(["step", "level"], [3, 10], [6, 11]), Select("date", ["2000-01-01"])) diff --git a/tests/test_slicing_unsliceable_axis.py b/tests/test_slicing_unsliceable_axis.py index 286252065..f5a73aab0 100644 --- a/tests/test_slicing_unsliceable_axis.py +++ b/tests/test_slicing_unsliceable_axis.py @@ -3,7 +3,6 @@ import pytest import xarray as xr -from polytope_feature.engine.hullslicer import HullSlicer from polytope_feature.polytope import Polytope, Request from polytope_feature.shapes import Box, Select from polytope_feature.utility.exceptions import UnsliceableShapeError @@ -17,9 +16,8 @@ def setup_method(self, method): dims=("date", "variable", "level"), coords={"date": pd.date_range("2000-01-01", "2000-01-03", 3), "variable": ["a"], "level": range(1, 130)}, ) - self.slicer = HullSlicer() options = {"compressed_axes_config": ["date", "variable", "level"]} - self.API = Polytope(datacube=array, engine=self.slicer, options=options) + self.API = Polytope(datacube=array, options=options) # Testing different shapes diff --git a/tests/test_slicing_xarray_3D.py b/tests/test_slicing_xarray_3D.py index 72f68d0ff..b1740810e 100644 --- a/tests/test_slicing_xarray_3D.py +++ b/tests/test_slicing_xarray_3D.py @@ -7,7 +7,6 @@ from polytope_feature.datacube.backends.xarray import XArrayDatacube from polytope_feature.datacube.tensor_index_tree import TensorIndexTree -from polytope_feature.engine.hullslicer import HullSlicer from polytope_feature.polytope import Polytope, Request from polytope_feature.shapes import ( Box, @@ -34,9 +33,8 @@ def setup_method(self, method): }, ) self.xarraydatacube = XArrayDatacube(array) - self.slicer = HullSlicer() options = {"compressed_axes_config": ["date", "step", "level"]} - self.API = Polytope(datacube=array, engine=self.slicer, options=options) + self.API = Polytope(datacube=array, options=options) # Testing different shapes diff --git a/tests/test_slicing_xarray_4D.py b/tests/test_slicing_xarray_4D.py index 9345e8db5..67db9292d 100644 --- a/tests/test_slicing_xarray_4D.py +++ b/tests/test_slicing_xarray_4D.py @@ -4,7 +4,6 @@ import xarray as xr from polytope_feature.datacube.tensor_index_tree import TensorIndexTree -from polytope_feature.engine.hullslicer import HullSlicer from polytope_feature.polytope import Polytope, Request from polytope_feature.shapes import ( Box, @@ -36,9 +35,8 @@ def setup_method(self, method): "lat": np.around(np.arange(0.0, 10.0, 0.1), 15), }, ) - self.slicer = HullSlicer() options = {"compressed_axes_config": ["date", "step", "level", "lat"]} - self.API = Polytope(datacube=array, engine=self.slicer, options=options) + self.API = Polytope(datacube=array, options=options) # Testing different shapes diff --git a/tests/test_snapping.py b/tests/test_snapping.py index 529f1d69a..240ef9729 100644 --- a/tests/test_snapping.py +++ b/tests/test_snapping.py @@ -1,7 +1,6 @@ import numpy as np import xarray as xr -from polytope_feature.engine.hullslicer import HullSlicer from polytope_feature.polytope import Polytope, Request from polytope_feature.shapes import Select @@ -17,9 +16,8 @@ def setup_method(self, method): "step": [1, 3, 5], }, ) - self.slicer = HullSlicer() options = {"compressed_axes_config": ["level", "step"]} - self.API = Polytope(datacube=array, engine=self.slicer, options=options) + self.API = Polytope(datacube=array, options=options) # Testing different shapes diff --git a/tests/test_snapping_real_data.py b/tests/test_snapping_real_data.py index 9b5950742..c2431d5c0 100644 --- a/tests/test_snapping_real_data.py +++ b/tests/test_snapping_real_data.py @@ -6,7 +6,6 @@ from earthkit import data from helper_functions import download_test_data -from polytope_feature.engine.hullslicer import HullSlicer from polytope_feature.polytope import Polytope, Request from polytope_feature.shapes import Box, Select @@ -17,20 +16,14 @@ def setup_method(self, method): download_test_data(nexus_url, "era5-levels-members.grib") ds = data.from_source("file", "./tests/data/era5-levels-members.grib") - array = ds.to_xarray(engine="cfgrib").isel(step=0).t - self.slicer = HullSlicer() - options = { + self.array = ds.to_xarray(engine="cfgrib").isel(step=0).t + self.options = { "axis_config": [ {"axis_name": "latitude", "transformations": [{"name": "reverse", "is_reverse": True}]}, {"axis_name": "longitude", "transformations": [{"name": "cyclic", "range": [0, 360]}]}, ], "compressed_axes_config": ["longitude", "latitude", "step", "time", "number", "isobaricInhPa"], } - self.API = Polytope( - datacube=array, - engine=self.slicer, - options=options, - ) @pytest.mark.internet def test_surrounding_on_grid_point(self): @@ -43,6 +36,11 @@ def test_surrounding_on_grid_point(self): Select("longitude", [requested_lon], method="surrounding"), Select("step", [np.timedelta64(0, "s")]), ) + + self.API = Polytope( + datacube=self.array, + options=self.options, + ) result = self.API.retrieve(request) result.pprint() country_points_plotting = [] diff --git a/tests/test_tree_protobuf_encoding.py b/tests/test_tree_protobuf_encoding.py index b177cbe2b..710c01694 100644 --- a/tests/test_tree_protobuf_encoding.py +++ b/tests/test_tree_protobuf_encoding.py @@ -49,17 +49,14 @@ def setup_method(self): def test_encoding(self): import pygribjump as gj - from polytope_feature.engine.hullslicer import HullSlicer from polytope_feature.polytope import Polytope self.options = { "pre_path": {"class": "od", "expver": "0001", "levtype": "sfc", "stream": "oper"}, } self.fdbdatacube = gj.GribJump() - self.slicer = HullSlicer() self.API = Polytope( datacube=self.fdbdatacube, - engine=self.slicer, options=self.options, ) fdb_datacube = self.API.datacube diff --git a/tests/test_tree_protobuf_encoding_fdb.py b/tests/test_tree_protobuf_encoding_fdb.py index 2e288f0d6..89a9075a2 100644 --- a/tests/test_tree_protobuf_encoding_fdb.py +++ b/tests/test_tree_protobuf_encoding_fdb.py @@ -12,7 +12,6 @@ def setup_method(self): def test_encoding(self): import pygribjump as gj - from polytope_feature.engine.hullslicer import HullSlicer from polytope_feature.polytope import Polytope, Request from polytope_feature.shapes import Box, Select @@ -61,10 +60,8 @@ def test_encoding(self): "pre_path": {"class": "od", "expver": "0001", "levtype": "sfc", "stream": "oper"}, } self.fdbdatacube = gj.GribJump() - self.slicer = HullSlicer() self.API = Polytope( datacube=self.fdbdatacube, - engine=self.slicer, options=self.options, ) result = self.API.retrieve(request) diff --git a/tests/test_type_change_transformation.py b/tests/test_type_change_transformation.py index 09f800610..2dd033f0d 100644 --- a/tests/test_type_change_transformation.py +++ b/tests/test_type_change_transformation.py @@ -1,7 +1,6 @@ import numpy as np import xarray as xr -from polytope_feature.engine.hullslicer import HullSlicer from polytope_feature.polytope import Polytope, Request from polytope_feature.shapes import Select @@ -21,8 +20,7 @@ def setup_method(self, method): "axis_config": [{"axis_name": "step", "transformations": [{"name": "type_change", "type": "int"}]}], "compressed_axes_config": ["step"], } - self.slicer = HullSlicer() - self.API = Polytope(datacube=array, engine=self.slicer, options=options) + self.API = Polytope(datacube=array, options=options) def test_merge_axis(self): request = Request(Select("step", [0])) diff --git a/tests/test_union_gj.py b/tests/test_union_gj.py index 9182fc360..13a99552d 100644 --- a/tests/test_union_gj.py +++ b/tests/test_union_gj.py @@ -1,7 +1,6 @@ import pandas as pd import pytest -from polytope_feature.engine.hullslicer import HullSlicer from polytope_feature.polytope import Polytope, Request from polytope_feature.shapes import Box, Select, Span, Union @@ -67,10 +66,8 @@ def test_fdb_datacube(self): ) self.fdbdatacube = gj.GribJump() - self.slicer = HullSlicer() self.API = Polytope( datacube=self.fdbdatacube, - engine=self.slicer, options=self.options, ) result = self.API.retrieve(request) @@ -108,10 +105,8 @@ def test_fdb_datacube_complete_overlap(self): ) self.fdbdatacube = gj.GribJump() - self.slicer = HullSlicer() self.API = Polytope( datacube=self.fdbdatacube, - engine=self.slicer, options=self.options, ) result = self.API.retrieve(request) @@ -149,10 +144,8 @@ def test_fdb_datacube_complete_overlap_v2(self): ) self.fdbdatacube = gj.GribJump() - self.slicer = HullSlicer() self.API = Polytope( datacube=self.fdbdatacube, - engine=self.slicer, options=self.options, ) result = self.API.retrieve(request) diff --git a/tests/test_union_point_box.py b/tests/test_union_point_box.py index ba280bc02..9132f14dc 100644 --- a/tests/test_union_point_box.py +++ b/tests/test_union_point_box.py @@ -1,7 +1,6 @@ import pandas as pd import pytest -from polytope_feature.engine.hullslicer import HullSlicer from polytope_feature.polytope import Polytope, Request from polytope_feature.shapes import Box, Point, Select, Union @@ -64,9 +63,7 @@ def test_fdb_datacube(self): Union(["latitude", "longitude"], box, point), ) self.fdbdatacube = gj.GribJump() - self.slicer = HullSlicer() self.API = Polytope( - request=request, datacube=self.fdbdatacube, engine=self.slicer, options=self.options,