diff --git a/NOTES b/NOTES index 36aaabc..4502eed 100644 --- a/NOTES +++ b/NOTES @@ -1,4 +1,40 @@ +* 2019-09-21 Version 0.8.0 + +Many improvements to pystare. +- Vertices & Centroids +- ConvexHull +- SpatialRange interval lists +- Expand intervals into spatial id values +- Intersection +- Visualization examples + +Some use is made of the STARE spatial index as a way to convey location information, +however there is some difficulty going back and forth to SpatialVectors. Note +that vertices and edges are not the easiest points to determine. + +There was a bug in expandIntervals due to the confusing nature of +EmbeddedLevelNameEncoding with its native and SciDB formats. + +Added a number of tests using Constraint. + +Improved intersection behavior by injecting dependency on varlen, i.e. +interior vs. boundary calculation. Setting varlen to false tends to +force calculation all the way down to leaf, whereas varlen true +stops on reaching a fully enclosed trixel, leading to the multi-resolution +partitioning. + +Corrected bug in SetPolyCorner of ConvexHull calculation. 2-corner list would +throw out co-linear third point without checking if it extended the side +or not. + +Changed HtmRangeMultiLevel::RangeFromIntersection to eliminate iteration over +the range in favor of using the SkipList's findMAX. Also eliminated superfluous +iteration after intersection is found. + +Added SpatialRange to use SkipList to generate intersections. Note the intersect +doesn't use find or search, so it's not clear we're getting the benefit of the SkipLists. + * 2019-09-05 Version 0.7.0 Added Griessenbaum's new apps and python interface. diff --git a/VERSION b/VERSION index bcaffe1..8adc70f 100644 --- a/VERSION +++ b/VERSION @@ -1 +1 @@ -0.7.0 \ No newline at end of file +0.8.0 \ No newline at end of file diff --git a/contrib/python/.gitignore b/contrib/python/.gitignore new file mode 100644 index 0000000..facc968 --- /dev/null +++ b/contrib/python/.gitignore @@ -0,0 +1 @@ +example-track-1.py diff --git a/contrib/python/example-1.py b/contrib/python/example-1.py new file mode 100644 index 0000000..463d8c1 --- /dev/null +++ b/contrib/python/example-1.py @@ -0,0 +1,8 @@ + +import numpy as np +import pystare as ps +a = np.array([0x0000000000000008, 0x000030000000000a, 0x000067ffffffffff, 0x000070000000000a, 0x0000907fffffffff],dtype=np.int64) +starts,ends = ps.from_intervals(a) +for i in range(starts.size): + print(i,hex(starts[i]),hex(ends[i])) + diff --git a/contrib/python/example-viz-1.py b/contrib/python/example-viz-1.py new file mode 100644 index 0000000..cebe7eb --- /dev/null +++ b/contrib/python/example-viz-1.py @@ -0,0 +1,217 @@ + +from math import ceil +import pystare as ps + +import matplotlib as mpl +import matplotlib.pyplot as plt +import matplotlib.tri as tri +import cartopy.crs as ccrs + +import numpy as np + +def shiftarg_lon(lon): + "If lon is outside +/-180, then correct back." + if(lon>180): + return ((lon + 180.0) % 360.0)-180.0 + else: + return lon + +def triangulate(i0,i1,i2): + "Prepare data structures for tri.Triangulate." + print('triangulating...') + # i0,i1,i2,ic = ps.to_vertices(indices) + i0lat,i0lon = ps.to_latlon(i0) + i1lat,i1lon = ps.to_latlon(i1) + i2lat,i2lon = ps.to_latlon(i2) + lats = np.zeros([3*len(i0lat)],dtype=np.double) + lons = np.zeros([3*len(i0lat)],dtype=np.double) + intmat = [] + k=0 + for i in range(len(i0)): + lats[k] = i0lat[i] + lons[k] = i0lon[i] + lats[k+1] = i1lat[i] + lons[k+1] = i1lon[i] + lats[k+2] = i2lat[i] + lons[k+2] = i2lon[i] + intmat.append([k,k+1,k+2]) + k=k+3 + for i in range(len(lons)): + lons[i] = shiftarg_lon(lons[i]) + print('triangulating done.') + return lons,lats,intmat + +def triangulate1(lats,lons): + "Prepare data for tri.Triangulate." + print('triangulating1...') + intmat=[] + npts=int(len(lats)/3) + k=0 + for i in range(npts): + intmat.append([k,k+1,k+2]) + k=k+3 + for i in range(len(lons)): + lons[i] = shiftarg_lon(lons[i]) + print('triangulating1 done.') + return lons,lats,intmat + +# Make a triangle's lat & lon corners +# lat = np.array([0, 0,60], dtype=np.double) +# lon = np.array([0,60,30], dtype=np.double) +lat = np.array([0, 0,60], dtype=np.double) +lon = np.array([0,60,30], dtype=np.double) + +# Make the indices out of the corners +s_resolution = 27 +indices = ps.from_latlon(lat, lon, s_resolution) +print('lat len: ',len(lat)) +print('indices len: ',len(indices)) +print('indices: ',[hex(i) for i in indices]) + +# Prepare to triangulate the indices (as locations) +def test1(indices): + "Plot the outline of a triangle using first three indices." + i0 = np.array([indices[0]], dtype=np.int64) + i1 = np.array([indices[1]], dtype=np.int64) + i2 = np.array([indices[2]], dtype=np.int64) + lons,lats,intmat = triangulate(i0,i1,i2) + print('test1 intmat len: ',len(intmat)) + print('test1 lats: ',lats) + print('test1 lons: ',lons) + print('test1 intm: ',intmat) + return lons,lats,intmat + # triang = tri.Triangulation(lons,lats,intmat) + # print('test1 triang : ',triang) + # return triang + +# Set up the projection and transformation +# proj = ccrs.PlateCarree() +# proj = ccrs.Robinson() +# proj = ccrs.Geodetic() +# proj = ccrs.Geodesic() +proj = ccrs.Mollweide() +transf = ccrs.Geodetic() +# transf = ccrs.PlateCarree() + +plt.figure() +plt.subplot(projection=proj,transform=transf) +# ax = plt.axes(projection=ccrs.PlateCarree()) +# ax = plt.axes(projection=ccrs.Mollweide()) +# proj=ccrs.Mollweide() +ax = plt.axes(projection=proj,transform=transf) +ax.set_global() +# ax.set_xlim(-180,180) +# ax.set_ylim(-90,90) +# ax.set_xlim(-1,1) +# ax.set_ylim(-1,1) +ax.coastlines() +# plt.contourf(xg,yg,v0g,60,transform=ccrs.PlateCarree()) +# plt.scatter(xg_flat,yg_flat,s=300,c=v_flat) +# plt.triplot(triang,'ko-') +# plt.show() + +def plot1(triang): + # plt.triplot(triang,'ro-',transform=ccrs.Geodetic()) + plt.triplot(triang,'r-',transform=transf) + # plt.show() + return + +test1_lons,test1_lats,test1_intmat = test1(indices) +plot1(tri.Triangulation(test1_lons,test1_lats,test1_intmat)) +plt.scatter(test1_lons,test1_lats,s=5,c='r',transform=transf) + +# resolution = 2; ntri=20 +resolution = 4; ntri=100 +# resolution = 7; ntri=1000 +# hull = ps.to_hull_range(indices,resolution,100) +hull = ps.to_hull_range_from_latlon(lat,lon,resolution,ntri) +print('0 hull len: ',len(hull)) + +# print(90) +# lats0 = np.zeros(len(hull)*4,dtype=np.int64) +# lons0 = np.zeros(len(hull)*4,dtype=np.int64) +# lats0,lons0 = ps._to_vertices_latlon(hull) +# print('lats0,lons0: ',len(lats0),len(lons0)) +# print(100) +lath,lonh,lathc,lonhc = ps.to_vertices_latlon(hull) +print('lath,lathc: ',len(lath),len(lathc)) +# print(110) +lons1,lats1,intmat1 = triangulate1(lath,lonh) + +## h0,h1,h2,hc = ps.to_vertices(hull) +## lons1,lats1,intmat1 = triangulate(h0,h1,h2) + +print('0 hull len: ',len(hull)) +print('0 hull lats len: ',len(lats1)) +jtest=18 +j=jtest +print('0 hull lats1: ',j,lats1[j*3:(j+1)*3]) +print('0 hull lons1: ',j,lons1[j*3:(j+1)*3]) +j=0 +# print('0 hull lats1: ',[i for i in lats1[j:j+12]]) +# print('0 hull lons1: ',[i for i in lons1[j:j+12]]) +print('') + +# tid = np.array([0x4c0000000000003],dtype=np.int64) +# t0,t1,t2,tc = ps.to_vertices(tid) +# print('t0: ',hex(t0[0])) +# print('t1: ',hex(t1[0])) +# print('t2: ',hex(t2[0])) +# print('tc: ',hex(tc[0])) +# print('t0 ll : ',ps.to_latlon(t0)) +# print('t1 ll : ',ps.to_latlon(t1)) +# print('t2 ll : ',ps.to_latlon(t2)) +# print('tc ll : ',ps.to_latlon(tc)) +# print('') + +if False: + # i=9; ilen=2 + # i=10; ilen=1 + # i=jtest; ilen=4 + i=jtest; ilen=10 + id_test = np.array(hull[i:i+ilen],dtype=np.int64) + print('i,id : ',i,[hex(j) for j in id_test]) + i0=i*3; i1=(i+ilen)*3 + lats1 = lats1[i0:i1] + lons1 = lons1[i0:i1] + # intmat1 = [intmat1[i]] + intmat1 = [] + for j in range(ilen): + intmat1.append([3*j,3*j+1,3*j+2]) + # intmat1 = [[0,2,1]] + id_test_lats,id_test_lons,id_test_latsc,id_test_lonsc = ps.to_vertices_latlon(id_test) + lonstest,latstest,intmattest = triangulate1(id_test_lats,id_test_lons) + # i0test,i1test,i2test,ictest = ps.to_vertices(id_test) + # print('test id: ',[hex(i) for i in id_test]) + # print('test ll : ',[(lats1[i],lons1[i]) for i in range(len(lats1))]) + # print('test im : ',[i for i in intmat1]) + # lonstest,latstest,intmattest = triangulate(i0test,i1test,i2test) + print('test lat: ',len(latstest)) + print('test lon: ',len(lonstest)) + print('test im: ',len(intmattest)) + print('test im[0]: ',intmattest[0]) + print('test im: ',intmattest) + triangtest = tri.Triangulation(lonstest,latstest,intmattest) + plt.triplot(triangtest,'g-',transform=transf) + plt.scatter(lonstest,latstest,s=5,c='g',transform=transf) + # plt.triplot(triangtest,'go-',transform=ccrs.Geodetic()) + # plt.scatter(lonstest,latstest,s=10,c='g',transform=ccrs.Geodetic()) + +print('resolution: ',resolution) +# print('hull : ',[hex(i) for i in hull]) +print('hull len: ',len(hull)) +print('hull ll len: ',len(lons1)) +print('hull im len: ',len(intmat1)) +# print('hull ll : ',[(lats1[i],lons1[i]) for i in range(len(lats1))]) +# print('hull im : ',[i for i in intmat1]) + +triang1 = tri.Triangulation(lons1,lats1,intmat1) +plt.triplot(triang1,'b-',transform=transf,lw=1,markersize=3) +plt.scatter(lons1,lats1,s=10,c='b',transform=transf) +# # plt.contourf(xg,yg,v0g,60,transform=ccrs.PlateCarree()) + +plt.show() + + + + diff --git a/contrib/python/example-viz-2.py b/contrib/python/example-viz-2.py new file mode 100644 index 0000000..53e3264 --- /dev/null +++ b/contrib/python/example-viz-2.py @@ -0,0 +1,98 @@ + +from math import ceil +import pystare as ps + +import matplotlib as mpl +import matplotlib.pyplot as plt +import matplotlib.tri as tri +import cartopy.crs as ccrs + +import numpy as np + +def shiftarg_lon(lon): + "If lon is outside +/-180, then correct back." + if(lon>180): + return ((lon + 180.0) % 360.0)-180.0 + else: + return lon + +def triangulate1(lats,lons): + "Prepare data for tri.Triangulate." + print('triangulating1...') + intmat=[] + npts=int(len(lats)/3) + k=0 + for i in range(npts): + intmat.append([k,k+1,k+2]) + k=k+3 + for i in range(len(lons)): + lons[i] = shiftarg_lon(lons[i]) + print('triangulating1 done.') + return lons,lats,intmat + +def plot1(lon,lat,lons,lats,triang,c0='r',c1='b',transf=None,lw=1): + if(lon is not None): + x=np.zeros([lon.size+1],dtype=np.double);x[:-1]=lon[:];x[-1]=lon[0] + y=np.zeros([lat.size+1],dtype=np.double); y[:-1]=lat[:]; y[-1]=lat[0] + ax.plot(x,y,True,transform=transf,c=c0) + plt.triplot(triang,c1+'-',transform=transf,lw=lw,markersize=3) + plt.scatter(lons,lats,s=10,c=c1,transform=transf) + return + +def make_hull(lat0,lon0,resolution0,ntri0): + hull0 = ps.to_hull_range_from_latlon(lat0,lon0,resolution0,ntri0) + lath0,lonh0,lathc0,lonhc0 = ps.to_vertices_latlon(hull0) + lons0,lats0,intmat0 = triangulate1(lath0,lonh0) + triang0 = tri.Triangulation(lons0,lats0,intmat0) + return lats0,lons0,triang0,hull0 + +resolution = 5 + +resolution0 = resolution; ntri0 = 1000 +lat0 = np.array([ 10, 5, 60,70], dtype=np.double) +lon0 = np.array([-30,-20,60,10], dtype=np.double) +lats0,lons0,triang0,hull0 = make_hull(lat0,lon0,resolution0,ntri0) +print('hull0: ',len(hull0)) + +resolution1 = resolution; ntri1 = 1000 +lat1 = np.array([10, 20, 30, 20 ], dtype=np.double) +lon1 = np.array([-60, 60, 60, -60], dtype=np.double) +lats1,lons1,triang1,hull1 = make_hull(lat1,lon1,resolution1,ntri1) +print('hull1: ',len(hull1)) + +if True: + intersected = np.full([1000],-1,dtype=np.int64) + intersected = ps.intersect(hull0,hull1,multiresolution=False) + # intersected = ps.intersect(hull0,hull1,multiresolution=True) + # print('hull0: ',[hex(i) for i in hull0]) + # print('hull1: ',[hex(i) for i in hull1]) + # ps._intersect_multiresolution(hull0,hull1,intersected) + # print('intersected: ',len(intersected)) + # print('np.min: ',np.amin(intersected)) + # print('intersected: ',[hex(i) for i in intersected]) + # The following are for _intersect_multiresolution's results + # endarg = np.argmax(intersected < 0) + # intersected = intersected[:endarg] + # intersected = ps.intersect(hull0,hull1) + print('intersected: ',len(intersected)) + lati,loni,latci,lonci = ps.to_vertices_latlon(intersected) + lonsi,latsi,intmati = triangulate1(lati,loni) + triangi = tri.Triangulation(lonsi,latsi,intmati) + +# Set up the projection and transformation +# proj = ccrs.PlateCarree() +# proj = ccrs.Robinson() +# proj = ccrs.Geodesic() +proj = ccrs.Mollweide() +transf = ccrs.Geodetic() +# transf = ccrs.PlateCarree() +plt.figure() +plt.subplot(projection=proj,transform=transf) +ax = plt.axes(projection=proj,transform=transf) +ax.set_global() +ax.coastlines() + +plot1(lon0,lat0,lons0,lats0,triang0,c0='r',c1='b',transf=transf) +plot1(lon1,lat1,lons1,lats1,triang1,c0='g',c1='c',transf=transf) +plot1(None,None,lonsi,latsi,triangi,c0='y',c1='r',transf=transf,lw=4) +plt.show() diff --git a/contrib/python/test.py b/contrib/python/test.py new file mode 100644 index 0000000..10b884a --- /dev/null +++ b/contrib/python/test.py @@ -0,0 +1,260 @@ +#!/usr/bin/python3 + +import numpy +import pystare + +lat = numpy.array([30,45,60], dtype=numpy.double) +lon = numpy.array([45,60,10], dtype=numpy.double) + +indices = pystare.from_latlon(lat, lon, 12) +print('0 indices: ',[hex(i) for i in indices]) + +lat, lon = pystare.to_latlon(indices) +print(lat, lon) + +lat, lon, level = pystare.to_latlonlevel(indices) +print(lat, lon, level) + +level = pystare.to_level(indices) +print(level) + +area = pystare.to_area(indices) +print(area) + +datetime = numpy.array(['2002-02-03T13:56:03.172', '2016-01-05T17:26:00.172'], dtype=numpy.datetime64) +print(datetime) + +index = pystare.from_utc(datetime, 6) +print(index) + +intersected = numpy.zeros([3], dtype=numpy.int64) +#? leni = 0 +pystare._intersect(indices, indices, intersected) +print('1 indices: ',[hex(i) for i in indices]) +print('1 intersected: ',[hex(i) for i in intersected] ) + +intersected = pystare.intersect(indices, indices) +print('2 intersected: ',[hex(i) for i in intersected]) + +intersected = pystare.intersect(indices, indices, multiresolution=True) +print('3 intersected: ',[hex(i) for i in intersected]) + +indices1 = numpy.array([indices[1]], dtype=numpy.int64) +intersected = pystare.intersect(indices, indices1, multiresolution=True) +print('4 intersected: ',[hex(i) for i in intersected]) + +indices1 = numpy.array([0x100000000000000c], dtype=numpy.int64) +intersected = pystare.intersect(indices, indices1, multiresolution=True) +print('5 intersected: ',[hex(i) for i in intersected]) +print('') + +indices1 = numpy.array([0,0,0], dtype=numpy.int64) +pystare._to_compressed_range(indices,indices1) +print('_indices1: ',[hex(i) for i in indices1]) + +indices1 = numpy.array([0,0,0], dtype=numpy.int64) +indices1 = pystare.to_compressed_range(indices) +print('indices1: ',list(map(hex,indices1))) + +indices1 = numpy.zeros([1000], dtype=numpy.int64) +result_size = numpy.zeros([1], dtype=numpy.int) +pystare._to_hull_range(indices,8,indices1,result_size) +# print('hull indices1: ',list(map(hex,indices1[0:100]))) +print('hull result_size: ',result_size) +# print('hull indices1: ',list(map(hex,indices1[899:910]))) +indices1=indices1[0:result_size[0]] +print('hull trimmed size ',indices1.size) +print('') + +hull_indices = pystare.to_hull_range(indices,8,2000) +print('hull_indices size: ',hull_indices.size) +print('') + +cmp = numpy.zeros([9],dtype=numpy.int64) +pystare._cmp_spatial(indices,indices,cmp) +print('cmp input: ',indices) +print('cmp: ',cmp,' i.e. diagonal') +print('') +cmp = numpy.zeros([3*1],dtype=numpy.int64) +indices1 = numpy.zeros([1],dtype=numpy.int64) +indices1[0] = indices[1] +print('cmp input1: ',indices) +print('cmp input2: ',indices1) +pystare._cmp_spatial(indices,indices1,cmp) +print('cmp: ',cmp) +print('') +cmp = numpy.zeros([3*2],dtype=numpy.int64) +indices1 = numpy.zeros([2],dtype=numpy.int64) +indices1[0] = indices[1] +indices1[1] = indices[1]+3 +print('cmp input1: ',indices) +print('cmp input2: ',indices1) +pystare._cmp_spatial(indices,indices1,cmp) +print('cmp: ',cmp) +pystare._cmp_spatial(indices1,indices,cmp) +print('cmp: ',cmp) +print('') +mp = numpy.zeros([3*2],dtype=numpy.int64) +indices1 = numpy.zeros([2],dtype=numpy.int64) +indices1[0] = indices[1]-2 +indices1[1] = indices[1] +print('cmp input1: ',indices) +print('cmp input2: ',indices1) +# pystare._cmp_spatial(indices,indices1,cmp) +cmp1 = pystare.cmp_spatial(indices,indices1) +print('cmp1: ',cmp1) +print('') + + +print('Expand intervals') +intervals_src = [\ + 0x2320000000000005,\ + 0x2324000000000005,\ + 0x2327ffffffffffff,\ + 0x3aa0000000000005,\ + 0x3aa7ffffffffffff,\ + 0x3aa8000000000004,\ + 0x3ab2000000000005,\ + 0x3ac7ffffffffffff,\ + 0x3ad0000000000005,\ + 0x3ae7ffffffffffff,\ + 0x3af2000000000005,\ + 0x3afa000000000005,\ + 0x3b40000000000004,\ + 0x3b4c000000000005,\ + 0x3b52000000000005,\ + 0x3b5a000000000005,\ + 0x3b5fffffffffffff,\ + 0x3b80000000000004,\ + 0x3b8a000000000005,\ + 0x3b8fffffffffffff,\ + 0x3b90000000000004,\ + 0x3b9fffffffffffff,\ + 0x3bc0000000000005,\ + 0x3bc7ffffffffffff,\ + 0x3bc8000000000004,\ + 0x3bd0000000000005,\ + 0x3bdfffffffffffff,\ + 0x3be2000000000005,\ + 0x3be8000000000004,\ + 0x3bf4000000000005,\ + 0x3bf8000000000005,\ + 0x3bfc000000000005,\ + 0x3bffffffffffffff,\ + 0x3e40000000000005,\ + 0x3e43ffffffffffff,\ + 0x3e46000000000005,\ + 0x3f20000000000005,\ + 0x3f24000000000005,\ + 0x3f27ffffffffffff,\ + 0x3f32000000000005,\ + 0x3f36000000000005,\ + 0x3f3a000000000005,\ + 0x3fa0000000000005\ +] + +intervals_expanded_src = [\ + 0x2320000000000005 ,\ + 0x2324000000000005 ,\ + 0x2326000000000005 ,\ + 0x3aa0000000000005 ,\ + 0x3aa2000000000005 ,\ + 0x3aa4000000000005 ,\ + 0x3aa6000000000005 ,\ + 0x3aa8000000000004 ,\ + 0x3ab2000000000005 ,\ + 0x3ab4000000000005 ,\ + 0x3ab6000000000005 ,\ + 0x3ab8000000000005 ,\ + 0x3aba000000000005 ,\ + 0x3abc000000000005 ,\ + 0x3abe000000000005 ,\ + 0x3ac0000000000005 ,\ + 0x3ac2000000000005 ,\ + 0x3ac4000000000005 ,\ + 0x3ac6000000000005 ,\ + 0x3ad0000000000005 ,\ + 0x3ad2000000000005 ,\ + 0x3ad4000000000005 ,\ + 0x3ad6000000000005 ,\ + 0x3ad8000000000005 ,\ + 0x3ada000000000005 ,\ + 0x3adc000000000005 ,\ + 0x3ade000000000005 ,\ + 0x3ae0000000000005 ,\ + 0x3ae2000000000005 ,\ + 0x3ae4000000000005 ,\ + 0x3ae6000000000005 ,\ + 0x3af2000000000005 ,\ + 0x3afa000000000005 ,\ + 0x3b40000000000004 ,\ + 0x3b4c000000000005 ,\ + 0x3b52000000000005 ,\ + 0x3b5a000000000005 ,\ + 0x3b5c000000000005 ,\ + 0x3b5e000000000005 ,\ + 0x3b80000000000004 ,\ + 0x3b8a000000000005 ,\ + 0x3b8c000000000005 ,\ + 0x3b8e000000000005 ,\ + 0x3b90000000000004 ,\ + 0x3b98000000000004 ,\ + 0x3bc0000000000005 ,\ + 0x3bc2000000000005 ,\ + 0x3bc4000000000005 ,\ + 0x3bc6000000000005 ,\ + 0x3bc8000000000004 ,\ + 0x3bd0000000000005 ,\ + 0x3bd2000000000005 ,\ + 0x3bd4000000000005 ,\ + 0x3bd6000000000005 ,\ + 0x3bd8000000000005 ,\ + 0x3bda000000000005 ,\ + 0x3bdc000000000005 ,\ + 0x3bde000000000005 ,\ + 0x3be2000000000005 ,\ + 0x3be8000000000004 ,\ + 0x3bf4000000000005 ,\ + 0x3bf8000000000005 ,\ + 0x3bfc000000000005 ,\ + 0x3bfe000000000005 ,\ + 0x3e40000000000005 ,\ + 0x3e42000000000005 ,\ + 0x3e46000000000005 ,\ + 0x3f20000000000005 ,\ + 0x3f24000000000005 ,\ + 0x3f26000000000005 ,\ + 0x3f32000000000005 ,\ + 0x3f36000000000005 ,\ + 0x3f3a000000000005 ,\ + 0x3fa0000000000005 \ +] +intervals = numpy.array(intervals_src,dtype=numpy.int64) +expected_expanded = numpy.array(intervals_expanded_src,dtype=numpy.int64) +intervals_expanded = numpy.zeros([len(expected_expanded)],dtype=numpy.int64) +expected_len = numpy.zeros([1],dtype=numpy.int64) + +print('len(intervals) ',len(intervals)) +print('len(expected_expanded) ',len(expected_expanded)) +print('len(intervals_expanded) ',len(intervals_expanded)) +print('len(expected_len) ',len(expected_len)) +intervals_len = len(intervals) +resolution = -1 +pystare._expand_intervals(intervals,resolution,intervals_expanded,expected_len) +# print('intervals_expanded: ',[hex(i) for i in intervals_expanded]) +print('expected_len: ',expected_len) +error_found = False +for i in range(len(intervals_expanded)): + if(intervals_expanded[i] != expected_expanded[i]): + error_found = True + print('_expanded_intervals error at i = ',i,' ',intervals_expanded[i],' != ',expected_expanded[i]) +if(not error_found): + print('_expanded_intervals ok') + +intervals_expanded_1 = pystare.expand_intervals(intervals,resolution) +for i in range(len(intervals_expanded_1)): + if(intervals_expanded_1[i] != expected_expanded[i]): + error_found = True + print('expanded_intervals error at i = ',i,' ',intervals_expanded_1[i],' != ',expected_expanded[i]) +if(not error_found): + print('expanded_intervals ok') diff --git a/include/CMakeLists.txt b/include/CMakeLists.txt index 352c1f1..b1a5c45 100644 --- a/include/CMakeLists.txt +++ b/include/CMakeLists.txt @@ -33,9 +33,10 @@ set ( ${CMAKE_CURRENT_SOURCE_DIR}/SpatialIndex.hxx ${CMAKE_CURRENT_SOURCE_DIR}/SpatialInterface.h ${CMAKE_CURRENT_SOURCE_DIR}/SpatialInterface.hxx - ${CMAKE_CURRENT_SOURCE_DIR}/SpatialSign.h + ${CMAKE_CURRENT_SOURCE_DIR}/SpatialSign.h ${CMAKE_CURRENT_SOURCE_DIR}/SpatialVector.h ${CMAKE_CURRENT_SOURCE_DIR}/SpatialVector.hxx + ${CMAKE_CURRENT_SOURCE_DIR}/SpatialRange.h ${CMAKE_CURRENT_SOURCE_DIR}/SpatialRotation.h ${CMAKE_CURRENT_SOURCE_DIR}/STARE.h ${CMAKE_CURRENT_SOURCE_DIR}/TemporalGeneral.h @@ -85,6 +86,7 @@ install(FILES STARE.h SpatialIndex.hxx SpatialInterface.h SpatialInterface.hxx + SpatialRange.h SpatialSign.h SpatialVector.h SpatialVector.hxx diff --git a/include/EmbeddedLevelNameEncoding.h b/include/EmbeddedLevelNameEncoding.h index d5adb91..15c9e1c 100644 --- a/include/EmbeddedLevelNameEncoding.h +++ b/include/EmbeddedLevelNameEncoding.h @@ -63,8 +63,12 @@ class EmbeddedLevelNameEncoding: public NameEncoding { uint32 aLevel ) const; ///< Somewhat dangerous to use. uint64 idFromTerminatorAndLevel_NoDepthBit(uint64 terminator, uint32 level); ///< Also a little dangerous. - bool terminatorp(); - bool terminatorp(uint64 terminatorCandidate); + + bool terminatorp() const; + bool terminatorp(uint64 terminatorCandidate) const; + + bool SciDBterminatorp() const; + bool SciDBterminatorp(uint64 terminatorCandidate) const; /// What triangle is just after the terminator? uint64 successorToTerminator_NoDepthBit(uint64 terminator, uint32 level) const; @@ -98,14 +102,15 @@ class EmbeddedLevelNameEncoding: public NameEncoding { EmbeddedLevelNameEncoding clearDeeperThanLevel(uint64 level); + void SciDBincrement_LevelToMaskDelta(uint32 level,uint64 &one_mask_to_level,uint64 &one_at_level) const; + void increment_LevelToMaskDelta(uint32 level,uint64 &one_mask_to_level,uint64 &one_at_level) const; + + /* + * increment and decrement work on the internal bit format, not the SciDB format. + */ uint64 increment(uint64 lowerBound, uint32 level, int steps = 1) const; uint64 decrement(uint64 lowerBound, uint32 level, int steps = 1) const; - bool terminatorp(uint64 terminator) const { - uint64 levelBits = terminator & levelMask; - return levelBits == levelMask; - } - EmbeddedLevelNameEncoding& operator=(EmbeddedLevelNameEncoding obj) { this->setId(obj.getId()); return *this; @@ -222,4 +227,6 @@ class EmbeddedLevelNameEncoding: public NameEncoding { }; +void EmbeddedLevelNameEncoding_test(); + #endif /* INCLUDE_EMBEDDEDLEVELNAMEENCODING_H_ */ diff --git a/include/HstmRange.h b/include/HstmRange.h index ec55f72..8d6db3e 100644 --- a/include/HstmRange.h +++ b/include/HstmRange.h @@ -36,6 +36,8 @@ namespace std { class HstmRange { public: HstmRange(); + HstmRange(HtmRangeMultiLevel_NameSpace::HtmRangeMultiLevel *range); + // TOO Deep copy? HstmRange(HstmRange range); virtual ~HstmRange(); // TODO Note the int in the following is a return code, not an index. @@ -53,8 +55,6 @@ class HstmRange { /** * For STARE we may need to coarsen position information to the resolution level since it points to a 7cm triangle inside a larger triangle. * - * - * */ void addRange(Key a, Key b); void addRange(HstmRange* r); diff --git a/include/HtmRangeMultiLevel.h b/include/HtmRangeMultiLevel.h index b3f0072..ab0f787 100644 --- a/include/HtmRangeMultiLevel.h +++ b/include/HtmRangeMultiLevel.h @@ -22,11 +22,11 @@ namespace HtmRangeMultiLevel_NameSpace { enum InclusionType { - InclOutside = 0, /* */ - InclInside, /* */ - InclLo, /* number is on low end of an interval */ - InclHi, /* number is on high end of an interval */ - InclAdjacentXXX + InclOutside = 0, /* */ + InclInside, /* */ + InclLo, /* number is on low end of an interval */ + InclHi, /* number is on high end of an interval */ + InclAdjacentXXX }; struct TInsideResult { @@ -47,89 +47,97 @@ struct TInsideResult { */ class LINKAGE HtmRangeMultiLevel { - public: - static int HIGHS; - static int LOWS; - static int BOTH; - int getNext(Key &lo, Key &hi); - int getNext(Key *lo, Key *hi); - KeyPair getNext(); - int getNext(KeyPair &kp); - - void setSymbolic(bool flag); - - void addRange(const Key lohi); - void addRange(const Key lo, const Key hi); - void addRange(HtmRangeMultiLevel *range); - void mergeRange(const Key lo, const Key hi); - void defrag(); - void defrag(Key gap); - void CompressionPass(); - void purge(); - int isIn(Key key); - int isIn(Key lo, Key hi); - int isIn(HtmRangeMultiLevel & otherRange); - Key bestgap(Key desiredSize); - - HtmRangeMultiLevel getSpan(); - void parse(std::string rangeString); - - bool equalp(HtmRangeMultiLevel *other); - - int stats(int desiredSize); - int nranges(); - int nindexes_in_ranges(); - - int getLosLength() { - return my_los->getLength(); - } - int getHisLength() { - return my_his->getLength(); - } - void reset(); - - void print(std::ostream& os, bool symbolic = false); // FIX THIS, so caller does not set symbolic here.... - void print(int what, std::ostream& os, bool symbolic = false); // FIX THIS, so caller does not set symbolic here.... - - int compare(const HtmRangeMultiLevel & other) const; - - /** Hold the encoding scheme for the symbolic representation. - */ - EmbeddedLevelNameEncoding *encoding; - /** Change or set symbolic encoding. - */ - void setEncoding(EmbeddedLevelNameEncoding *encoding) { this->encoding = encoding; } - /** Return the encoding for translations, etc. - */ - EmbeddedLevelNameEncoding *getEncoding() { return encoding; } - - // Moved here per ajmendez. - // TODO MLR Why two skip lists for ranges? Why not one skip list? Or an interval arithmetic package? - SkipList *my_los; - SkipList *my_his; - bool symbolicOutput; - - HtmRangeMultiLevel(); - HtmRangeMultiLevel(EmbeddedLevelNameEncoding *encoding); - ~HtmRangeMultiLevel(){ - purge(); - delete encoding; - delete my_los; - delete my_his; - }; - - HtmRangeMultiLevel *HtmRangeMultiLevelAtLevelFromIntersection(HtmRangeMultiLevel *range2, int htmIdLevel=-1); - int contains(Key a, Key b); - - friend LINKAGE ostream& operator<<(ostream& os, const HtmRangeMultiLevel& range); - - protected: - TInsideResult tinside(const Key mid) const; - // const char buffer[256]; -// bool symbolicOutput; -// private: -// SkipList *my_los; -// SkipList *my_his; +public: + + // For the stored indices + bool embeddedLevel = true; // Set to false for deprecated legacy behavior. + + // For the range data type + static int HIGHS; + static int LOWS; + static int BOTH; + int getNext(Key &lo, Key &hi); + int getNext(Key *lo, Key *hi); + KeyPair getNext(); + int getNext(KeyPair &kp); + + void setSymbolic(bool flag); + + void addRange(const Key lohi); + void addRange(const Key lo, const Key hi); + void addRange(HtmRangeMultiLevel *range); + void mergeRange(const Key lo, const Key hi); + void defrag(); + void defrag(Key gap); + void CompressionPass(); + void purge(); + int isIn(Key key); + int isIn(Key lo, Key hi); + int isIn(HtmRangeMultiLevel & otherRange); + Key bestgap(Key desiredSize); + + HtmRangeMultiLevel getSpan(); + void parse(std::string rangeString); + + bool equalp(HtmRangeMultiLevel *other); + + int stats(int desiredSize); + int nranges(); + int nindexes_in_ranges(); + + int getLosLength() { + return my_los->getLength(); + } + int getHisLength() { + return my_his->getLength(); + } + void reset(); + + void print(std::ostream& os, bool symbolic = false); // FIX THIS, so caller does not set symbolic here.... + void print(int what, std::ostream& os, bool symbolic = false); // FIX THIS, so caller does not set symbolic here.... + + int compare(const HtmRangeMultiLevel & other) const; + + /** Hold the encoding scheme for the symbolic representation. + */ + EmbeddedLevelNameEncoding *encoding; + /** Change or set symbolic encoding. + */ + void setEncoding(EmbeddedLevelNameEncoding *encoding) { this->encoding = encoding; } + /** Return the encoding for translations, etc. + */ + EmbeddedLevelNameEncoding *getEncoding() { return encoding; } + + // Moved here per ajmendez. + // TODO MLR Why two skip lists for ranges? Why not one skip list? Or an interval arithmetic package? + SkipList *my_los; + SkipList *my_his; + bool symbolicOutput; + + HtmRangeMultiLevel(); + HtmRangeMultiLevel(EmbeddedLevelNameEncoding *encoding); + ~HtmRangeMultiLevel(){ + purge(); + delete encoding; + delete my_los; + delete my_his; + }; + + HtmRangeMultiLevel *HtmRangeMultiLevelAtLevelFromIntersection(HtmRangeMultiLevel *range2, int htmIdLevel=-1); + HtmRangeMultiLevel *RangeFromIntersection(HtmRangeMultiLevel *range2, bool compress = true, int force_htmIdLevel=-1); + + int contains(Key a, Key b); + // TODO KeyPair contains(Key a, Key b); + + friend LINKAGE ostream& operator<<(ostream& os, const HtmRangeMultiLevel& range); + +protected: + TInsideResult tinside(const Key mid) const; + // const char buffer[256]; + // bool symbolicOutput; + // private: + // SkipList *my_los; + // SkipList *my_his; }; #define SKIP_PROB (0.5) diff --git a/include/PySTARE.h b/include/PySTARE.h index 133a3b4..a5876c4 100644 --- a/include/PySTARE.h +++ b/include/PySTARE.h @@ -25,17 +25,31 @@ #include #include "STARE.h" +#include "SpatialRange.h" static STARE stare; // Spatial -void from_latlon(double* lat, int len_lat, double * lon, int len_lon, int64_t* indices, int level); +void from_latlon(double* lat, int len_lat, double * lon, int len_lon, int64_t* indices, int level); void to_latlon(int64_t* indices, int len, double* lat, double* lon); void to_latlonlevel(int64_t* indices, int len, double* lat, double* lon, int* levels); void to_level(int64_t* indices, int len, int* levels); void to_triangle(int64_t* indices, int len); +// Broken or dangerous. void to_vertices(int64_t* indices, int len, int64_t* vertices0, int64_t* vertices1, int64_t* vertices2, int64_t* centroid); +void _to_vertices_latlon(int64_t* indices, int len, double* triangle_info_lats, int dmy1, double* triangle_info_lons, int dmy2 ); void to_area(int64_t* indices, int len, double* areas); +void _to_compressed_range(int64_t* indices, int len, int64_t* range_indices, int len_ri); +void _to_hull_range (int64_t* indices, int len, int resolution, int64_t* range_indices, int len_ri, int64_t* result_size, int len_rs); +void _expand_intervals(int64_t* indices, int len, int resolution, int64_t* range_indices, int len_ri, int64_t* result_size, int len_rs); + +void _to_hull_range_from_latlon(double* lat, int len_lat, double* lon, int len_lon, int resolution, int64_t* range_indices, int len_ri, int64_t* result_size, int len_rs); +void from_intervals(int64_t* intervals, int len, int64_t* indices_starts, int64_t* indices_terminators ); + +void _intersect(int64_t* indices1, int len1, int64_t* indices2, int len2, int64_t* intersection, int leni); +void _intersect_multiresolution(int64_t* indices1, int len1, int64_t* indices2, int len2, int64_t* intersection, int leni); +void _cmp_spatial(int64_t* indices1, int len1, int64_t* indices2, int len2, int64_t* cmp, int len12); + // Temporal void from_utc(int64_t *datetime, int len, int64_t *indices, int resolution); //void to_utc(int64_t* indices, int len, double* julian_day); diff --git a/include/STARE.h.in b/include/STARE.h.in index 1be7eaa..a39e252 100644 --- a/include/STARE.h.in +++ b/include/STARE.h.in @@ -43,8 +43,11 @@ typedef std::vector STARE_ArrayIndexSpatialValues; /// TODO The following is a stopgap on the way to completely removing this hardcoded bit. #define STARE_HARDWIRED_RESOLUTION_LEVEL_MAX 27 +// < 2019-0919 #define STARE_HARDWIRED_BUILD_LEVEL_MAX 5 #define STARE_HARDWIRED_BUILD_LEVEL_MAX 5 - +// #define STARE_HARDWIRED_BUILD_LEVEL_MAX 7 +// #define STARE_HARDWIRED_BUILD_LEVEL_MAX 8 +// #define STARE_HARDWIRED_BUILD_LEVEL_MAX 10 bad alloc class STARE { public: @@ -60,12 +63,17 @@ public: // Spatial array index functions. [Maybe change the name StareId to spatialStareId in the following...] STARE_ArrayIndexSpatialValue ValueFromLatLonDegrees(float64 latDegrees, float64 lonDegrees, int resolutionLevel = STARE_HARDWIRED_RESOLUTION_LEVEL_MAX); LatLonDegrees64 LatLonDegreesFromValue(STARE_ArrayIndexSpatialValue spatialStareId); + SpatialVector SpatialVectorFromValue(STARE_ArrayIndexSpatialValue spatialStareId); + STARE_ArrayIndexSpatialValue ValueFromSpatialVector(SpatialVector v, int resolution); uint32 ResolutionLevelFromValue(STARE_ArrayIndexSpatialValue spatialStareId); Triangle TriangleFromValue(STARE_ArrayIndexSpatialValue spatialStareId, int resolutionLevel = -1); float64 AreaFromValue (STARE_ArrayIndexSpatialValue spatialStareId, int resolutionLevel = -1); + + STARE_ArrayIndexSpatialValues toVertices(STARE_ArrayIndexSpatialValues spatialStareIds); STARE_SpatialIntervals ConvexHull(LatLonDegrees64ValueVector points,int force_resolution_level); + STARE_SpatialIntervals ConvexHull(STARE_ArrayIndexSpatialValues points,int force_resolution_level); SpatialIndex getIndex() { return sIndex; } SpatialIndex getIndex(int resolutionLevel); @@ -143,7 +151,22 @@ STARE_ArrayIndexSpatialValue sTerminator(STARE_ArrayIndexSpatialValue spatialSta bool terminatorp(STARE_ArrayIndexSpatialValue spatialStareId); /// Check if the index value is a terminator. -HstmRange SpatialRangeFromSpatialIntervals(STARE_SpatialIntervals intervals); +// Deprecated... cf. SpatialRange - HstmRange SpatialRangeFromSpatialIntervals(STARE_SpatialIntervals intervals); + +void STARE_ArrayIndexSpatialValues_insert(STARE_ArrayIndexSpatialValues& v, STARE_ArrayIndexSpatialValue siv); + +STARE_SpatialIntervals expandInterval(STARE_SpatialIntervals interval, int64 force_resolution = -1); +STARE_ArrayIndexSpatialValues expandIntervals(STARE_SpatialIntervals intervals, int64 force_resolution = -1); + +STARE_ArrayIndexSpatialValue shiftSpatialIdAtLevel( + STARE_ArrayIndexSpatialValue spatialStareId, + int resolution, + int shiftAmount + ); + +STARE_SpatialIntervals spatialIntervalFromHtmIDKeyPair(KeyPair kp); + +uint64 spatialLevelMask(); void STARE_test(); diff --git a/include/SkipListElement.h b/include/SkipListElement.h index 062bec2..a072100 100644 --- a/include/SkipListElement.h +++ b/include/SkipListElement.h @@ -36,7 +36,8 @@ typedef int64 Key; // key type TODO Why not uint64? // typedef uint64 Key; // key type TODO Why not uint64? -typedef int Value; // value type +// typedef int Value; // value type +typedef int64 Value; // mlr... // class ostream; class SkipListElement; diff --git a/include/SpatialConstraint.h b/include/SpatialConstraint.h index de6e908..109fe9e 100644 --- a/include/SpatialConstraint.h +++ b/include/SpatialConstraint.h @@ -81,6 +81,10 @@ class LINKAGE SpatialConstraint : public SpatialSign { /// Initialization constructor SpatialConstraint(SpatialVector, float64); + /// Construct a zero constraint through two points + // Makes a zero-angle constraint directed towards v0^v1. + SpatialConstraint(SpatialVector v0, SpatialVector v1); + /// Copy constructor SpatialConstraint(const SpatialConstraint &); diff --git a/include/SpatialGeneral.h b/include/SpatialGeneral.h index fe8644d..1f952ab 100644 --- a/include/SpatialGeneral.h +++ b/include/SpatialGeneral.h @@ -115,6 +115,8 @@ const float64 gPi = 3.1415926535897932385E0 ; const float64 gPio2 = 3.1415926535897932385E0/2.0 ; const float64 gPr = gPi/180.0; const float64 piDiv180 = gPr; + +// TODO Set these by environment variable or other property list. // const float64 gEpsilon = 1.0E-15; // const float64 gEpsilon = 1.0E-16; // const float64 gEpsilon = 1.0E-17; diff --git a/include/SpatialInterface.h b/include/SpatialInterface.h index c0bd648..047a16b 100644 --- a/include/SpatialInterface.h +++ b/include/SpatialInterface.h @@ -48,10 +48,11 @@ struct htmPolyCorner { /** htmInterface class. + The SpatialInterface class contains all methods to interface the HTM index with external applications. - // TODO Consider renaming htmInterface to SpatialInterface + // TODO Consider renaming htmInterface to SpatialInterface, but maybe "htm" is a good indicator that we're working at a low level with this interface. */ @@ -72,6 +73,8 @@ class LINKAGE htmInterface { saveDepth parameter can be specified to keep the given amount of levels in memory. This can also be altered by changeDepth. */ + htmInterface(const SpatialIndex *index); + htmInterface(size_t searchlevel = 5, size_t buildevel = 5, SpatialRotation rot = rot_identity ); // [ed:gyuri:saveDepth was 2] /** Destructor. */ @@ -187,7 +190,8 @@ class LINKAGE htmInterface { * @return */ const HTMRangeValueVector & convexHull( LatLonDegrees64ValueVector latlon, - size_t steps = -1); + size_t steps = -1, + bool interiorp = false); /** Request all triangles in the convex hull of a given set of points. @@ -302,6 +306,7 @@ class LINKAGE htmInterface { // void setPolyCorner(SpatialVector &v); // this routine does the work for all convexHull calls + bool hull_interiorp_ = false; const HTMRangeValueVector & doHull(); }; diff --git a/include/SpatialRange.h b/include/SpatialRange.h new file mode 100644 index 0000000..485ee0d --- /dev/null +++ b/include/SpatialRange.h @@ -0,0 +1,57 @@ +/* + * SpatialRange.h + * + * Created on: Sep 8, 2019 + * Author: mrilee + * + * Copyright (C) 2019 Rilee Systems Technologies LLC + */ + +#ifndef INCLUDE_SPATIALRANGE_H_ +#define INCLUDE_SPATIALRANGE_H_ + +#include "STARE.h" +#include "HstmRange.h" + +/** + * A wrapper for an HstmRange that knows about STARE to provide htm-range-like functions. + * + */ +class SpatialRange { +public: + SpatialRange(); + SpatialRange(STARE_SpatialIntervals intervals); + SpatialRange(HstmRange *range) { this->range = range; } + virtual ~SpatialRange(); + + void addSpatialIntervals(STARE_SpatialIntervals intervals); + void addSpatialRange(SpatialRange r); + + STARE_SpatialIntervals toSpatialIntervals(); + int getNextSpatialInterval(STARE_SpatialIntervals &interval); + + HstmRange *range; + + /////////// private: Maybe? //////////// + + // Maybe inherit these? + // TODO Note the int in the following is a return code, not an index. + int getNext(KeyPair &kp) { + int istat = range->getNext(kp); +// cout << "" << flush; + return istat; + }; + void reset() { range->reset(); } // range not null? + void purge() { range->purge(); } // what if range null? +}; + +SpatialRange* sr_intersect(const SpatialRange& a, const SpatialRange& b, bool compress = false); + +inline SpatialRange* operator& ( const SpatialRange& a, const SpatialRange& b) { + return sr_intersect(a,b); + // return new SpatialRange(new HstmRange(a.range->range->RangeFromIntersection(b.range->range))); // NOTE mlr Probably about the safest way to inst. SpatialRange. +} +void SpatialRange_test(); + + +#endif /* INCLUDE_SPATIALRANGE_H_ */ diff --git a/src/CMakeLists.txt b/src/CMakeLists.txt index f636b08..d5fa51c 100644 --- a/src/CMakeLists.txt +++ b/src/CMakeLists.txt @@ -26,6 +26,7 @@ set ( SpatialVector.cpp SpatialRotation.C STARE.C + SpatialRange.cpp TemporalIndex.cpp TemporalWordFormat1.cpp VarStr.C @@ -79,9 +80,9 @@ if (SWIG_FOUND AND POLICY CMP0078 AND POLICY CMP0086) # swig_add_library( PySTARE TYPE SHARED LANGUAGE python SOURCES PySTARE.i ) swig_add_library( pystare TYPE SHARED LANGUAGE python SOURCES PySTARE.i ) - swig_link_libraries( pystare STARE ) swig_link_libraries( pystare PySTAREA ) swig_link_libraries( pystare ${PYTHON_LIBRARIES} ) + swig_link_libraries( pystare STARE ) endif() install(TARGETS STARE DESTINATION lib) diff --git a/src/EmbeddedLevelNameEncoding.C b/src/EmbeddedLevelNameEncoding.C index 5d97ce3..0fdf1f4 100644 --- a/src/EmbeddedLevelNameEncoding.C +++ b/src/EmbeddedLevelNameEncoding.C @@ -9,6 +9,7 @@ #include #include #include +#include EmbeddedLevelNameEncoding::EmbeddedLevelNameEncoding() {} @@ -22,11 +23,13 @@ EmbeddedLevelNameEncoding::~EmbeddedLevelNameEncoding() {} * @return */ char* EmbeddedLevelNameEncoding::nameById(uint64 id) { +#ifdef DIAG if(id == 0) { // Throw an exception? std::cout << "EmbeddedLevelNameEncoding::nameById WARNING ID == 0 -> 'S0'..."; // throw SpatialFailure("EmbeddedLevelNameEncoding::nameById-INVALID_ID_0"); } +#endif int nameSize = levelById(id)+3; ///< levelById is local to the encoding char *returnedName = new char[nameSize]; @@ -190,7 +193,7 @@ uint64 EmbeddedLevelNameEncoding::idFromTerminatorAndLevel_NoDepthBit(uint64 ter int64 EmbeddedLevelNameEncoding::getSciDBLeftJustifiedFormat(uint64 leftId) const { uint64 id_NoLevelBit = leftId & stripLevelBitMask; // Note this covers bits 0-5. - uint64 level = leftId & levelMask; + uint64 level = leftId & levelMask; // TODO Repent the sin of redundant code. int64 leftId_scidb = id_NoLevelBit; @@ -221,7 +224,7 @@ int64 EmbeddedLevelNameEncoding::getSciDBLeftJustifiedFormat() const { // *** TODO *** NEED ERROR CHECKING & SIGNALLING FOR INVALID IDs - int64 leftId_scidb = 0; // Erf... // What's the invalid id now? Anything negative. + int64 leftId_scidb = 0; // Erf... // What's the invalid id now? Anything negative? uint64 id_NoLevelBit = this->maskOffLevelBit(); int64 level = this->getLevel(); @@ -238,6 +241,9 @@ int64 EmbeddedLevelNameEncoding::getSciDBLeftJustifiedFormat() const { } int64 EmbeddedLevelNameEncoding::getSciDBTerminatorLeftJustifiedFormat() const { + if(terminatorp()) { + return getSciDBLeftJustifiedFormat(); // Just return what we have. + } return getSciDBLeftJustifiedFormat(getIdTerminator_NoDepthBit()); } @@ -248,8 +254,12 @@ void EmbeddedLevelNameEncoding::setIdFromSciDBLeftJustifiedFormat( int64 id_scid uint64 level = id_scidb & levelMaskSciDB; iTmp = iTmp << 1; iTmp = iTmp | level; - iTmp = iTmp | TopBit; + iTmp = iTmp | TopBit; // This says we have a valid ID. How does this relate to HTM-ID? // this->id = iTmp; // hacking... + // if(SciDBterminatorp(id_scidb)) { + if( level == levelMaskSciDB ) { + iTmp = iTmp | levelMask; // Fill in the extra bit. + } this->setId(iTmp); } @@ -353,20 +363,52 @@ uint64 EmbeddedLevelNameEncoding::predecessorToLowerBound_NoDepthBit(uint64 lowe return terminator; } -bool EmbeddedLevelNameEncoding::terminatorp() { +bool EmbeddedLevelNameEncoding::terminatorp() const { return terminatorp(this->id); } +bool EmbeddedLevelNameEncoding::terminatorp(uint64 terminator) const { + uint64 levelBits = terminator & levelMask; + return levelBits == levelMask; +} +bool EmbeddedLevelNameEncoding::SciDBterminatorp() const { + return SciDBterminatorp(this->id); +} +bool EmbeddedLevelNameEncoding::SciDBterminatorp(uint64 terminator) const { + uint64 levelBits = terminator & levelMaskSciDB; + return levelBits == levelMaskSciDB; +} + +void EmbeddedLevelNameEncoding::increment_LevelToMaskDelta(uint32 level,uint64 &one_mask_to_level,uint64 &one_at_level) const { + one_mask_to_level = 0; + one_at_level = one; + for(uint64 shift = 2; + shift <= ( (topBitPosition-3) - 2*(level) ); + shift +=2 ) { + one_mask_to_level = one_mask_to_level << 2; + one_mask_to_level += 3; + one_at_level = one_at_level << 2; + } +// one_at_level = one_at_level >> 2; +} + +void EmbeddedLevelNameEncoding::SciDBincrement_LevelToMaskDelta(uint32 level,uint64 &one_mask_to_level,uint64 &one_at_level) const { + increment_LevelToMaskDelta(level,one_mask_to_level,one_at_level); + one_mask_to_level = one_mask_to_level >> 1; + one_at_level = one_at_level >> 1; -bool EmbeddedLevelNameEncoding::terminatorp(uint64 terminatorCandidate) { - uint64 level = terminatorCandidate & levelMask; - return level == 63; +// one_at_level = one_at_level >> 2; } uint64 EmbeddedLevelNameEncoding::increment(uint64 lowerBound, uint32 level, int n) const { /// TODO Error checking of overflow not trustworthy here. using namespace std; + + bool topBitSet = (lowerBound & TopBit) == TopBit; // Some left-justified may not have a topbit set. + uint64 successor = lowerBound; // Bump up one, but we still need the level. // Should clean up successor just in case terminator non-3 prefix is not consistent with level. + + /* uint64 one_mask_to_level = 0; uint64 one_at_level = one; for(uint64 shift = 2; @@ -377,18 +419,62 @@ uint64 EmbeddedLevelNameEncoding::increment(uint64 lowerBound, uint32 level, int one_at_level = one_at_level << 2; } // one_at_level = one_at_level >> 2; +*/ - successor = successor & (~one_mask_to_level); + uint64 one_mask_to_level, one_at_level; + increment_LevelToMaskDelta(level,one_mask_to_level,one_at_level); + +#ifdef DIAG + cout << "lowerBound: " << setw(23) << hex << lowerBound << endl << flush; + cout << "successor: " << setw(23) << hex << successor << endl << flush; +#endif + + // Remove TopBit to test for overflow. + successor = successor & (~TopBit); + + // Truncate to level + successor = successor & (~one_mask_to_level); + + // Increment successor += n*one_at_level; +#ifdef DIAG + cout << "successor: " << setw(23) << hex << successor << endl << flush; + cout << "n: " << setw(23) << dec << n << endl << flush; + cout << "level: " << setw(23) << dec << level << endl << flush; + cout << "n*one_at_level: " << setw(23) << hex << n*one_at_level << endl << flush; + cout << "one_at_level: " << setw(23) << hex << one_at_level << endl << flush; + cout << "successor&TB: " << setw(23) << hex << (successor & TopBit) << endl << flush; + cout << "TopBit: " << setw(23) << hex << TopBit << endl << flush; +#endif + // cout << "one_at_level: "<< hex << one_at_level << dec << endl << flush; // cout << "one_mask_to_: "<< hex << one_mask_to_level << dec << endl << flush; - // Check for overflow. - if( successor == TopBit ) { - return 0; // It's invalid! + // Check for overflow. Not a completely accurate check. + if( (successor & TopBit) == TopBit ) { + +#ifdef DIAG + cout << "n: " << setw(23) << dec << n << endl << flush; + cout << "lowerBound: " << setw(23) << hex << lowerBound << endl << flush; + cout << "level: " << setw(23) << dec << level << endl << flush; + cout << "successor: " << setw(23) << hex << successor << endl << flush; + cout << "n*one_at_level: " << setw(23) << hex << n*one_at_level << endl << flush; + cout << "one_at_level: " << setw(23) << hex << one_at_level << endl << flush; + cout << "TopBit: " << setw(23) << hex << TopBit << endl << flush; +#endif + + throw SpatialFailure("EmbeddedLevelNameEncoding::error-increment-overflow"); + // TODO THROW EXCEPTION INSTEAD! + // cout << "ERROR 1" << endl << flush; + // cout << "successor " << successor << ", topbit " << TopBit << endl << flush; + // return 0; // It's invalid! } + // Add TopBit back in if necessary. + if( topBitSet ) { + successor = successor | TopBit; + } successor += level; @@ -403,9 +489,11 @@ uint64 EmbeddedLevelNameEncoding::increment(uint64 lowerBound, uint32 level, int // if( successor < (( lowerBound & stripMask ) + level) ) { // if( (successor & ~levelMask) < (( lowerBound & ~levelMask )) ) { if( (successor & stripMask) < ( lowerBound & stripMask ) ) { - return 0; // It's invalid! Wrap around! + throw SpatialFailure("EmbeddedLevelNameEncoding::error-increment-wrap-around"); + // TODO THROW EXCEPTION INSTEAD! + // cout << "ERROR 2" << endl << flush; + // return 0; // It's invalid! Wrap around! } - return successor; } uint64 EmbeddedLevelNameEncoding::decrement(uint64 lowerBound, uint32 level, int n) const { @@ -432,7 +520,9 @@ uint64 EmbeddedLevelNameEncoding::decrement(uint64 lowerBound, uint32 level, int successor -= n*one_at_level; if( successor == 0 ) { - return 0; // It's invalid! + throw SpatialFailure("EmbeddedLevelNameEncoding::error-decrement-underflow"); + // TODO THROW EXCEPTION INSTEAD! + // return 0; // It's invalid! } successor += level; @@ -453,10 +543,14 @@ uint64 EmbeddedLevelNameEncoding::decrement(uint64 lowerBound, uint32 level, int // Check for underflow if( (successor & stripMask) > (lowerBound & stripMask)) { + throw SpatialFailure("EmbeddedLevelNameEncoding::error-decrement-wrap-around"); // if( successor > (lowerBound & ~levelMask)+level) { - return 0; // It's invalid! Wrap around! + // TODO THROW EXCEPTION INSTEAD! + // return 0; // It's invalid! Wrap around! } return successor; } + + diff --git a/src/HstmRange.C b/src/HstmRange.C index 19ff0ba..4f82c9b 100644 --- a/src/HstmRange.C +++ b/src/HstmRange.C @@ -6,16 +6,19 @@ */ #include "HstmRange.h" +#include namespace std { HstmRange::HstmRange() { - // TODO Auto-generated constructor stub range = new HtmRangeMultiLevel_NameSpace::HtmRangeMultiLevel(); } +HstmRange::HstmRange(HtmRangeMultiLevel_NameSpace::HtmRangeMultiLevel *range){ + this->range = range; +} + HstmRange::~HstmRange() { - // TODO Auto-generated destructor stub delete range; } @@ -25,7 +28,7 @@ HstmRange::~HstmRange() { * @param b_ a Key (int64) */ void HstmRange::addRange(Key a_, Key b_) { - Key a = leftJustifiedEncoding.maskOffLevelBit(a_); + Key a = leftJustifiedEncoding.maskOffLevelBit(a_); // The level bit is a bit at the top used as a validity check. A deprecated feature and not needed for a left-justified index value. int aLevel = leftJustifiedEncoding.levelById(a_); Key b = leftJustifiedEncoding.maskOffLevelAndLevelBit(b_); int bLevel = leftJustifiedEncoding.levelById(b_); @@ -35,13 +38,22 @@ void HstmRange::addRange(Key a_, Key b_) { // cout << dec << 902 << " : aLevel = " << aLevel << endl << flush; if( aLevel != bLevel ) { - cout << "HstmRange::addRange::ERROR::KeyLevelMismatch " - << aLevel << " != " << bLevel << " " - << "a = " << a_ << " " - << "b = " << b_ << " " - << endl << flush; - // TODO Throw? - throw SpatialException("HstmRange::addRange::ERROR::KeyLevelMismatch"); + if( !leftJustifiedEncoding.terminatorp(b_) ) { +#ifdef DIAG + cout << "HstmRange::addRange::ERROR::KeyLevelMismatch " + << aLevel << " != " << bLevel << " " + << "a = " << a_ << " " + << "b = " << b_ << " " + << endl << flush; +#endif + // TODO Throw? + throw SpatialException("HstmRange::addRange::ERROR::KeyLevelMismatch"); + } else { + b = b | leftJustifiedEncoding.levelMask; + } + } else { + leftJustifiedEncoding.setId(b_); + b = leftJustifiedEncoding.getIdTerminator_NoDepthBit(); } // Put a into HtmRange without change @@ -65,8 +77,20 @@ void HstmRange::addRange(Key a_, Key b_) { Key bSave = b; */ + /* 2019-0909 MLR THIS PIECE "WORKS". leftJustifiedEncoding.setId(b_); b = leftJustifiedEncoding.getIdTerminator_NoDepthBit(); + */ + +// #define DIAG +#undef DIAG +#define IDOUT(p,m,s) p << m << " " << setw(16) << setfill('0') << hex << s << dec << " " << s << endl << flush; +#ifdef DIAG + IDOUT(cout,"hr::ar a_: ",a_) + IDOUT(cout,"hr::ar a : ",a) + IDOUT(cout,"hr::ar b_: ",b_) + IDOUT(cout,"hr::ar b : ",b) +#endif // cout << dec << 1000 << " : " << hex << "0x" << a_ << " 0x" << bSave << endl << flush; diff --git a/src/HtmRangeMultiLevel.cpp b/src/HtmRangeMultiLevel.cpp index b21cd40..a67ddfe 100644 --- a/src/HtmRangeMultiLevel.cpp +++ b/src/HtmRangeMultiLevel.cpp @@ -12,7 +12,7 @@ #include // setw() #include // various *RepresentationString elements #include -#include // levelOfId +#include // levelOfId // Do wwe really need this? // TODO If levelOfId not needed, remove. #ifdef _WIN32 #include @@ -27,6 +27,15 @@ using namespace std; using namespace HtmRangeMultiLevel_NameSpace; +uint64 HRML_levelOfId(uint64 id,bool embeddedLevel,uint64 levelMask) { + if( embeddedLevel ) { + return id & levelMask; + } else { + // legacy behavior + return levelOfId(id); + } +} + /** * Translate an HtmRangeMultiLevel to one at a greater level. If the desired level is * less that the level implicit in the input range (lo & hi), then just return @@ -41,11 +50,11 @@ using namespace HtmRangeMultiLevel_NameSpace; */ KeyPair HtmRangeMultiLevelAtLevelFromHtmRangeMultiLevel(int htmIdLevel, Key lo, Key hi) { // htmIdLevel is used to set maxlevel in the index. aka olevel. - int levelLo = levelOfId(lo); + int levelLo = HRML_levelOfId(lo,false,63); if(levelLonranges()<=0)||(range2->nranges()<=0)) return 0; @@ -78,14 +87,14 @@ HtmRangeMultiLevel *HtmRangeMultiLevel::HtmRangeMultiLevelAtLevelFromIntersectio if (!indexp1) return 0; if(htmIdLevel<0) { - htmIdLevel = levelOfId(lo1); + htmIdLevel = HRML_levelOfId(lo1,false,63); } // cout << "indexp1: " << indexp1 << endl << flush; // cout << "l,lo,hi1: " << htmIdLevel << " " << lo1 << " " << hi1 << endl << flush; - // cout << "a" << flush; + cout << "a" << flush; do { - // cout << "b" << endl << flush; + cout << "b" << endl << flush; KeyPair testRange1 = HtmRangeMultiLevelAtLevelFromHtmRangeMultiLevel(htmIdLevel,lo1,hi1); range2->reset(); uint64 indexp2 = range2->getNext(lo2,hi2); @@ -119,8 +128,111 @@ HtmRangeMultiLevel *HtmRangeMultiLevel::HtmRangeMultiLevelAtLevelFromIntersectio return resultRange; } +KeyPair HRML_AtLevelFromMultiLevel(uint64 htmIdLevel, Key lo, Key hi, uint64 levelMask) { + uint64 levelLo = lo & levelMask; + const uint64 one = 1; + if( levelLo < htmIdLevel) { + lo = (lo & ~levelMask) | htmIdLevel; + // Ignore weird cases where level(lo) != level(hi) + // Make hi into a terminator for htmIdLevel. + hi = hi | ( (one << (6+54-2*htmIdLevel)) - 1 ); + } + KeyPair levelAdaptedRange; + levelAdaptedRange.lo = lo; + levelAdaptedRange.hi = hi; + return levelAdaptedRange; +} + +/** + * Find intersection of two ranges and return as a range. + * + * TODO Replace double-nested-loop with SkipLists search or find functionality. + * + */ +HtmRangeMultiLevel *HtmRangeMultiLevel::RangeFromIntersection( + HtmRangeMultiLevel *range2, bool compress, int force_htmIdLevel ) { + HtmRangeMultiLevel *range1 = this; // Just an alias + if((!range1)||(!range2)) return 0; + if((range1->nranges()<=0)||(range2->nranges()<=0)) return 0; + Key lo1,hi1,lo2,hi2; + range1->reset(); + uint64 indexp1 = range1->getNext(lo1,hi1); + if(!indexp1) return 0; + if(force_htmIdLevel<0) { + force_htmIdLevel = lo1 & this->encoding->levelMask; // TODO Establish 31 or 63? + } + HtmRangeMultiLevel *resultRange = new HtmRangeMultiLevel(); resultRange->purge(); + + do { + KeyPair testRange1 = HRML_AtLevelFromMultiLevel(force_htmIdLevel,lo1,hi1,this->encoding->levelMask); + range2->reset(); // Sigh. Reset and loop from the beginning. TODO Avoid restarting loop. There must be a faster way. + /* Try to skip past by using SkipList functions. */ + // TODO do something like range2->findMAX... + Key loKey = range2->my_los->findMAX(testRange1.lo); + Value hiKey = range2->my_los->search(loKey,true); + Value vhi = range2->my_his->search(hiKey,true); + /**/ + uint64 indexp2 = range2->getNext(lo2,hi2); // TODO Implement a find or search for inserting. + bool intersects = false, past_chance; +#define FMTX(x) setw(16) << setfill('0') << hex << x << dec + // Search forward until we find an intersection. Once an intersection is found, + // figure out what the intersection is and add it to the result range. + // int kount=0; + if(indexp2) + do { + // ++kount; + KeyPair testRange2 = HRML_AtLevelFromMultiLevel(force_htmIdLevel,lo2,hi2,this->encoding->levelMask); + intersects = testRange2.lo <= testRange1.hi + && testRange2.hi >= testRange1.lo; +// #define DIAG +#ifdef DIAG + cout << "lh1,lh2: " + << FMTX(lo1) << " " << FMTX(hi1) << ", " + << FMTX(lo2) << " " << FMTX(hi2) << ", " + << intersects << flush; +#endif + if(intersects){ + Key lo_ = max(testRange1.lo,testRange2.lo); + Key hi_ = min(testRange1.hi,testRange2.hi); + resultRange->addRange(lo_,hi_); +#ifdef DIAG + cout << ", added " + << FMTX(lo_) << " " + << FMTX(hi_) << flush; +#endif + } +#ifdef DIAG + cout << "." << endl << flush; +#endif +#undef DIAG + past_chance = (uint64) testRange2.lo > (uint64) testRange1.hi; + } while (range2->getNext(lo2,hi2) && !past_chance); + // cout << "kount = " << kount << endl << flush; + } while (range1->getNext(lo1,hi1)); // TODO Can we replace getNext with some sort of find or search. +#undef FMTX + // cout << "d" << flush; + // cout << "d nr " << resultRange->nranges() << endl << flush; + // cout << "d rr " << hex << resultRange << dec << endl << flush; + if(resultRange->nranges()>0) { + if(compress) { + resultRange->CompressionPass(); + } + resultRange->defrag(); + } + // cout << "e" << flush; + return resultRange; +} + + // Note the default use of EmbeddedLevelNameEncoding. -HtmRangeMultiLevel::HtmRangeMultiLevel() : HtmRangeMultiLevel(new EmbeddedLevelNameEncoding()) {} +// HtmRangeMultiLevel::HtmRangeMultiLevel() : HtmRangeMultiLevel(new EmbeddedLevelNameEncoding()) {} +HtmRangeMultiLevel::HtmRangeMultiLevel() { + encoding = new EmbeddedLevelNameEncoding(); + my_los = new SkipList(SKIP_PROB); + // cout << "hrml my_los " << hex << my_los << dec << endl << flush; + my_his = new SkipList(SKIP_PROB); + symbolicOutput = false; +} HtmRangeMultiLevel::HtmRangeMultiLevel(EmbeddedLevelNameEncoding *encoding) { this->encoding = encoding; @@ -449,7 +561,8 @@ void HtmRangeMultiLevel::mergeRange(const Key lo, const Key hi) // Add the first one. if( my_los->myHeader->getElement(0) == NIL ) { - my_los->insert(lo,100); + // my_los->insert(lo,100); + my_los->insert(lo,hi); my_his->insert(hi,100); // cout << 8002 << " First one inserted. " << endl << flush; return; @@ -508,7 +621,8 @@ void HtmRangeMultiLevel::mergeRange(const Key lo, const Key hi) // Don't know what's above h. Iterate. } else if( hi1 < l ) { // Case 1. A is below B. Just add - my_los->insert(lo1,10001); + // my_los->insert(lo1,10001); + my_los->insert(lo1,hi1); my_his->insert(hi1,10001); done = true; } else if( (lo1 < l) && ( (l <= hi1) && (hi1 <= h) ) ) { @@ -537,14 +651,16 @@ void HtmRangeMultiLevel::mergeRange(const Key lo, const Key hi) // At the same level, merge the two. my_los->freeRange(l_m,h_p); // Freeing up to h_p is okay because h_p==h is still part of current interval. my_his->freeRange(l_m,h_p); - my_los->insert(l_m,100021); + // my_los->insert(l_m,100021); + my_los->insert(l_m,h_p); my_his->insert(h_p,100021); } else { // The lower part overlaps an empty part. Just add. // ??? Don't need a freeRange ??? Okay... if(level > l_level) { // cout << "8000-1031" << endl << flush; // If the new interval's level is greater, just skip in the current, add before. - my_los->insert(l_m,100022); + // my_los->insert(l_m,100022); + my_los->insert(l_m,h_m); my_his->insert(h_m,100022); } else if(true) { // Case 2.3 level < l_level -- new interval wins // TODO WORRY -- What about collisions? If we have a collision, will we simply put in the value back in? @@ -561,10 +677,12 @@ void HtmRangeMultiLevel::mergeRange(const Key lo, const Key hi) // } my_los->freeRange(l_m,h_p); my_his->freeRange(l_m,h_p); - my_los->insert(l_m,100023); // lo1 // this changes the level to level + // my_los->insert(l_m,100023); // lo1 // this changes the level to level + my_los->insert(l_m,h_0); // lo1 // this changes the level to level my_his->insert(h_0,100023); // hi1 if(h_0 < h_p) { - my_los->insert(l_p,100023); + // my_los->insert(l_p,100023); + my_los->insert(l_p,h_p); my_his->insert(h_p,100023); } my_los->reset(); my_his->reset(); @@ -597,7 +715,8 @@ void HtmRangeMultiLevel::mergeRange(const Key lo, const Key hi) cout << "HtmRangeMultiLevel::mergeRange::ERROR!!! SUCC(H) > HI1, I.E. THEY'RE EQUIVALENT." << endl << flush; } if( l_m < h_m ) { // If lo1 and l are equivalent, current one wins, and we ignore the new one. - my_los->insert(l_m,1000421); + // my_los->insert(l_m,1000421); + my_los->insert(l_m,h_m); my_his->insert(h_m,1000421); // my_los->reset(); my_his->reset(); // } else { @@ -633,13 +752,16 @@ void HtmRangeMultiLevel::mergeRange(const Key lo, const Key hi) my_los->freeRange(l_m,h_p); my_his->freeRange(l_m,h_p); if( l_m < h_m ) { // If l_m and l_0 are "equivalent", so the current interval wins and we ignore interval_m. - my_los->insert(l_m,10004); + // my_los->insert(l_m,10004); + my_los->insert(l_m,h_m); my_his->insert(h_m,10004); } - my_los->insert(l_0,10004); + // my_los->insert(l_0,10004); + my_los->insert(l_0,h_0); my_his->insert(h_0,10004); // TODO Subtle bug? Need to verify edge case. if( l_p < h_p ) { // Non equivalent h_0 and h_p. - my_los->insert(l_p,10004); + // my_los->insert(l_p,10004); + my_los->insert(l_p,h_p); my_his->insert(h_p,10004); } } @@ -694,10 +816,12 @@ void HtmRangeMultiLevel::mergeRange(const Key lo, const Key hi) my_his->freeRange(l_m,h_0); // TODO NOTE: When predecessor is used, you have to check to see if pred(b) is less than inf(interval). if( l_m < h_m ) { // If l_m and l_0 are "equivalent", so the current interval wins and we ignore interval_m. - my_los->insert(l_m,100052); + // my_los->insert(l_m,100052); + my_los->insert(l_m,h_m); my_his->insert(h_m,100052); } - my_los->insert(l_0,100052); + // my_los->insert(l_0,100052); + my_los->insert(l_0,h_0); my_his->insert(h_0,100052); // TODO Subtle bug? Need to verify edge case. // update for iteration lo1 = l_p; @@ -715,7 +839,8 @@ void HtmRangeMultiLevel::mergeRange(const Key lo, const Key hi) if(not done) { // The new interval goes at the end of the skiplists. // Case 6. We're at the top. Just add. - my_los->insert(lo1,10006); + // my_los->insert(lo1,10006); + my_los->insert(lo1,hi1); my_his->insert(hi1,10006); done = true; } @@ -741,7 +866,7 @@ void HtmRangeMultiLevel::mergeRange(const Key lo, const Key hi) */ void HtmRangeMultiLevel::addRange(const Key lo, const Key hi) { -// my_los->insert(lo, (Value) 0); // TODO Consider doing something useful with (Value)... +// my_los->insert(lo, (Value) 0); // TODO Consider doing something useful with (Value)... Like storing hi... // my_his->insert(hi, (Value) 0); // cout << "x200: " << hex << lo << " " << hi << endl; // cout << "x201: " << (lo == hi) << endl; @@ -1010,28 +1135,57 @@ void HtmRangeMultiLevel::CompressionPass() { Key newLoPredecessor = encoding->predecessorToLowerBound_NoDepthBit(newLo,level0); // oldLo..newLoPredecessor; newLo..hi0 Modify the skiplists. my_his->insert(newLoPredecessor,1024); - my_los->insert(newLo,1024); + // my_los->insert(newLo,1024); + my_los->insert(newLo,newLoPredecessor); // Set lists to the new lo my_los->reset(); my_his->reset(); // TODO Until we know better, start over. Bad, bad, bad. // TODO Perhaps instead try a find or a search that would set the iterators. } else { // cout << "400: " << endl << flush; +// cout << "400: blo = 0x" << setw(16) << setfill('0') << hex << bareLo << endl << flush; +// cout << "400: bhi = 0x" << setw(16) << setfill('0') << hex << bareHi << endl << flush; +// cout << "400: dlt = 0x" << setw(16) << setfill('0') << hex << delta << endl << flush; +// cout << "400: dlt = " << setw(16) << dec << delta << endl << flush; // Snip off as much as possible int numberToCoalesce = (delta+1) / 4; + // cout << "410: ntc " << dec << numberToCoalesce << endl << flush; Key oldLo = lo0; Key newLo = oldLo; - for(int i=0; iincrement(newLo,level0); - } - } + // cout << "420: lo0 = 0x" << setw(16) << setfill('0') << hex << lo0 << endl << flush; + +// for(int i=0; iincrement(newLo,level0); +// } +// } + my_los->free(oldLo); --oldLo; // Reduce level - my_los->insert(oldLo,1025); - Key newLoPredecessor = encoding->predecessorToLowerBound_NoDepthBit(newLo,level0); - if(newLoPredecessor != hi0) { - my_his->insert(newLoPredecessor,1025); - my_los->insert(newLo,1025); + // my_los->insert(oldLo,1025); + my_los->insert(oldLo,my_his->getkey()); // What's the current key (for my_his)? // TODO Don't worry. + + try { + // Scoop up a bunch + for(int i=0; iincrement(newLo,level0); + } + } + // Then break it in half. + Key newLoPredecessor = encoding->predecessorToLowerBound_NoDepthBit(newLo,level0); + if(newLoPredecessor != hi0) { + my_his->insert(newLoPredecessor,1025); + // my_los->insert(newLo,1025); + my_los->insert(newLo,newLoPredecessor); + } + } catch ( SpatialException &e ) { + // What if we're at the top index already? Ooops, there's no new low at the new break. + // cout << "400: " << e.what() << endl << flush; + if( string(e.what()) != string("EmbeddedLevelNameEncoding::error-increment-overflow") ) { + throw SpatialFailure("HtmRangeMultiLevel::Compress::unknown error while incrementing to newLo."); + } } + + my_los->reset(); my_his->reset(); // TODO Reset is too drastic. Prefer to step back a little... Bad, bad, bad. } } @@ -1076,14 +1230,17 @@ void HtmRangeMultiLevel::reset() /// The number of ranges. int HtmRangeMultiLevel::nranges() { -// cout << "z000" << endl << flush; + // cout << "z000" << endl << flush; Key lo; // Key hi; int n_ranges; n_ranges = 0; + // cout << "z001" << endl << flush; + // cout << "z001 my_los " << hex << my_los << dec << endl << flush; my_los->reset(); + // cout << "z002" << endl << flush; my_his->reset(); -// cout << "z010" << endl << flush; + // cout << "z010" << endl << flush; // This is a problem when lo can be zero. Is it? // getkey returns -1 if nothing is found, maybe fix the following using >= 0? Worry about id 0. Should be okay this low in the code. MLR 2019-0327 @@ -1581,8 +1738,10 @@ int HtmRangeMultiLevel::getNext(Key &lo, Key &hi) // cout << " " << hi << " " << flush; // OLD if (hi <= (Key) 0){ if (hi < (Key) 0){ +#if DIAG cout << endl; cout << " getNext error!! " << endl << flush; +#endif // hi = lo = (Key) 0; hi = lo = (Key) -1; return 0; diff --git a/src/NOTES b/src/NOTES index 23b228e..23d2238 100644 --- a/src/NOTES +++ b/src/NOTES @@ -1,3 +1,14 @@ + + + +2019-0920 HtmRangeMultiLevel Compress bug + +TODO add a unit test for edge cases for this file +The compress bug resulted from the logic not accounting for what happens when you're at the last index value at a level. I.e. you're at the end of the road. + + +2019 Early + test if a triangle (id) is rejected: Implicitly, id is used in macro NV diff --git a/src/PySTARE.cpp b/src/PySTARE.cpp index 990fdc6..bc17f7d 100644 --- a/src/PySTARE.cpp +++ b/src/PySTARE.cpp @@ -37,8 +37,7 @@ void to_latlonlevel(int64_t* indices, int len, double* lat, double* lon, int* le void to_level(int64_t* indices, int len, int* levels) { for (int i=0; i0) { cout << "0x" << setw(16) << setfill('0') << hex << result[i] << " "; } + range_indices[i] = result[i]; + } + // cout << dec << endl << flush; + result_size[0] = result.size(); +} + +void _to_hull_range_from_latlon(double* lat, int len_lat, double* lon, int len_lon, int resolution, int64_t* range_indices, int len_ri, int64_t* result_size, int len_rs) { + + LatLonDegrees64ValueVector points; + for(int i=0; i0) { cout << "0x" << setw(16) << setfill('0') << hex << result[i] << " "; } + range_indices[i] = result[i]; + } + // cout << dec << endl << flush; + result_size[0] = result.size(); +} + +void _intersect(int64_t* indices1, int len1, int64_t* indices2, int len2, int64_t* intersection, int leni) { + STARE_SpatialIntervals si1(indices1, indices1+len1), si2(indices2, indices2+len2); + // intersection[0] = 69; + SpatialRange r1(si1), r2(si2); + // SpatialRange *ri = r1 & r2; + SpatialRange *ri = sr_intersect(r1,r2,false); + STARE_SpatialIntervals result_intervals = ri->toSpatialIntervals(); + STARE_ArrayIndexSpatialValues result = expandIntervals(result_intervals); + leni = result.size(); + for(int i=0; itoSpatialIntervals(); + STARE_ArrayIndexSpatialValues result = expandIntervals(result_intervals); + // cout << 400 << endl << flush; + leni = result.size(); + // cout << 500 << " result size " << result.size() << endl << flush; + for(int i=0; i 0): + endarg = endarg + 1 + # endarg = numpy.argmax(range_indices < 0) + range_indices = range_indices[:endarg] + return range_indices + +def expand_intervals(intervals,resolution,result_size_limit=1000): + result = numpy.full([result_size_limit],-1,dtype=numpy.int64) + result_size = numpy.full([1],-1,dtype=numpy.int64) + _expand_intervals(intervals,resolution,result,result_size) + result = result[:result_size[0]] + return result + +def to_hull_range(indices,resolution,range_size_limit=1000): + out_length = range_size_limit + range_indices = numpy.full([out_length],-1,dtype=numpy.int64) + result_size = numpy.full([1],-1,dtype=numpy.int64) + _to_hull_range(indices,resolution,range_indices,result_size) + range_indices = range_indices[:result_size[0]] + return range_indices + +def to_hull_range_from_latlon(lat,lon,resolution,range_size_limit=1000): + out_length = range_size_limit + range_indices = numpy.full([out_length],-1,dtype=numpy.int64) + result_size = numpy.full([1],-1,dtype=numpy.int64) + _to_hull_range_from_latlon(lat,lon,resolution,range_indices,result_size) + range_indices = range_indices[:result_size[0]] + return range_indices + +def to_vertices_latlon(indices): + out_length = len(indices) + lats = numpy.zeros([4*out_length],dtype=numpy.double) + lons = numpy.zeros([4*out_length],dtype=numpy.double) + # _to_vertices_latlon(indices,lats,lons,0) + lats,lons = _to_vertices_latlon(indices) + latsv = numpy.zeros([3*out_length],dtype=numpy.double) + lonsv = numpy.zeros([3*out_length],dtype=numpy.double) + latc = numpy.zeros([out_length],dtype=numpy.double) + lonc = numpy.zeros([out_length],dtype=numpy.double) + + k=0; l=0 + for i in range(out_length): + latsv[l] = lats[ k ] + lonsv[l] = lons[ k ] + + latsv[l+1] = lats[ k+1 ] + lonsv[l+1] = lons[ k+1 ] + + latsv[l+2] = lats[ k+2 ] + lonsv[l+2] = lons[ k+2 ] + + latc [i] = lats[ k+3 ] + lonc [i] = lons[ k+3 ] + k = k + 4; l = l + 3 + return latsv,lonsv,latc,lonc + +def cmp_spatial(indices1, indices2): + out_length = len(indices1)*len(indices2) + cmp = numpy.zeros([out_length],dtype=numpy.int64) + _cmp_spatial(indices1,indices2,cmp) + return cmp + +def intersect(indices1, indices2, multiresolution=True): + out_length = 2*max(len(indices1),len(indices2)) + intersected = numpy.full([out_length],-1,dtype=numpy.int64) + leni = 0 + if(multiresolution): + _intersect_multiresolution(indices1, indices2, intersected) + else: + _intersect(indices1, indices2, intersected) + endarg = numpy.argmax(intersected < 0) + intersected = intersected[:endarg] + return intersected +%} + %include "PySTARE.h" diff --git a/src/RangeConvex.cpp b/src/RangeConvex.cpp index b9915e8..db5f8a7 100644 --- a/src/RangeConvex.cpp +++ b/src/RangeConvex.cpp @@ -920,7 +920,7 @@ SpatialMarkup RangeConvex::testTrixel(uint64 nodeIndex) { // cout << endl << flush; // cout << " check children & partials " << endl << flush; - // NEW NEW algorithm Disabled when enablenew is 0 + // NEW < algorithm Disabled when enablenew is 0 { childID = indexNode->childID_[0]; if (childID != 0) { @@ -951,7 +951,7 @@ SpatialMarkup RangeConvex::testTrixel(uint64 nodeIndex) { If partial, then we look ahead to see how many children are rejected. But ah, next iteration could benefit from having computed this already. - If two chidlren are rejected, then we stop + If two children are rejected, then we stop If one or 0 nodes are rejected, then we */ // cout << " mark: " << mark << endl << flush; diff --git a/src/STARE.C b/src/STARE.C index 3a25943..2310424 100644 --- a/src/STARE.C +++ b/src/STARE.C @@ -14,6 +14,13 @@ #include "SpatialDomain.h" #include "SpatialInterface.h" #include +#include + +#ifndef DIAG +#define DIAGOUT1(a) +#else +#define DIAGOUT1(a) a +#endif /** * @brief Version function with C linkage to aid in finding the library with autoconf @@ -22,6 +29,21 @@ extern "C" const char *STARE_version() { return (const char *)STARE_VERSION; } + +void STARE_ArrayIndexSpatialValues_insert(STARE_ArrayIndexSpatialValues& v, STARE_ArrayIndexSpatialValue siv) { + // cout << dec << 1000 << endl; + STARE_ArrayIndexSpatialValues::iterator + it = std::lower_bound( v.begin(), v.end(), siv); + // cout << dec << 1010 << endl; + if(it == v.end()) { + v.insert( it, siv ); + } else if( (*it) != siv) { + // cout << dec << 1020 << endl; + v.insert( it, siv ); + // cout << dec << 1030 << endl; + } + // cout << dec << 1040 << endl; +} /** * @@ -114,6 +136,34 @@ STARE_ArrayIndexSpatialValue STARE::ValueFromLatLonDegrees( return leftJustifiedWithResolution.getSciDBLeftJustifiedFormat(); } +STARE_ArrayIndexSpatialValue STARE::ValueFromSpatialVector(SpatialVector v, int resolution) { + uint64 htmID; + SpatialVector vtry(v); + int k = 3; + while(k>0) { + try { + --k; + htmID = sIndex.idByPoint(vtry); + } catch( SpatialException e ) { + cerr << e.what(); + if( k > 0 ) { + cerr << " " << k << " Trying again... " << endl << flush; + vtry = vtry + 1.0e-10*SpatialVector(1,0,0); vtry.normalize(); + } else { + cerr << endl << flush; + stringstream ss; ss << setprecision(16); + ss << "STARE::ValueFromSpatialVector can't find vector v= " << v << endl; + throw SpatialFailure(ss.str().c_str()); + } + + } + } + BitShiftNameEncoding rightJustified(htmID); + EmbeddedLevelNameEncoding leftJustified(rightJustified.leftJustifiedId()); + EmbeddedLevelNameEncoding leftJustifiedWithResolution = leftJustified.atLevel(resolution, true); // True means keep all bits + return leftJustifiedWithResolution.getSciDBLeftJustifiedFormat(); +} + /** * Extract the resolution information from the spatial array index value. Since this * does not use the sIndex, it doesn't really need to be a method of this class. @@ -141,6 +191,7 @@ LatLonDegrees64 STARE::LatLonDegreesFromValue(STARE_ArrayIndexSpatialValue spati // cout << "sid: " << spatialStareId << endl << flush; uint64 htmID = htmIDFromValue(spatialStareId); + // cout << "lldfv htmID " << setw(16) << setfill('0') << hex << htmID << dec << endl << flush; SpatialVector v; /// This returns the center of the triangle (at index.search_level). Need to extract the position information. @@ -149,11 +200,13 @@ LatLonDegrees64 STARE::LatLonDegreesFromValue(STARE_ArrayIndexSpatialValue spati float64 lat=-999, lon=-999; v.getLatLonDegrees(lat, lon); - // cout << "sid-latlon: " << lat << ", " << lon << endl << flush; + // cout << "0 sid-latlon: " << lat << ", " << lon << endl << flush; // LatLonDegrees64 latlon = {.lat = lat, .lon = lon }; LatLonDegrees64 latlon(lat, lon); // = {.lat = lat, .lon = lon }; + // cout << "1 sid-latlon: " << latlon.lat << ", " << latlon.lon << endl << flush; + // return latlon; return LatLonDegrees64(lat, lon); @@ -165,6 +218,16 @@ LatLonDegrees64 STARE::LatLonDegreesFromValue(STARE_ArrayIndexSpatialValue spati } +SpatialVector STARE::SpatialVectorFromValue(STARE_ArrayIndexSpatialValue spatialStareId) { + // uint64 htmID = htmIDFromValue(spatialStareId,STARE_HARDWIRED_RESOLUTION_LEVEL_MAX); // Max resolution + uint64 htmID = htmIDFromValue(spatialStareId); // Max resolution + // cout << "svfv htmID " << setw(16) << setfill('0') << hex << htmID << dec << endl << flush; + SpatialVector v; + /// This returns the center of the triangle (at index.search_level). Need to extract the position information. + sIndex.pointByHtmId(v, htmID); + return v; +} + Triangle STARE::TriangleFromValue(STARE_ArrayIndexSpatialValue spatialStareId, int resolutionLevel) { // Users are going to expect the default resolution level to be that embedded in the sStareId. uint64 htmID = -1; @@ -212,6 +275,17 @@ Triangle STARE::TriangleFromValue(STARE_ArrayIndexSpatialValue spatialStareId, i return {.centroid=vc, .vertices=vertices}; } +STARE_ArrayIndexSpatialValues STARE::toVertices(STARE_ArrayIndexSpatialValues spatialStareIds) { + STARE_ArrayIndexSpatialValues spatialValues; + for(int i=0; i -1 ) { - cout << dec << 1100 << endl << flush; + // cout << dec << 1100 << endl << flush; // htm = htmInterface(&index_); htm = new htmInterface( this->getIndex(force_resolution_level).getMaxlevel(), this->getIndex(force_resolution_level).getBuildLevel(), this->getIndex(force_resolution_level).getRotation()); - cout << dec << 1101 << endl << flush; + // cout << dec << 1101 << endl << flush; } else { - cout << dec << 1200 << endl << flush; + // cout << dec << 1200 << endl << flush; // htm = htmInterface(&index_); htm = new htmInterface( this->getIndex(8).getMaxlevel(), this->getIndex(8).getBuildLevel(), this->getIndex(8).getRotation()); - cout << dec << 1201 << endl << flush; + // cout << dec << 1201 << endl << flush; } - cout << dec << "a2000" << endl << flush; + // cout << dec << "a2000" << endl << flush; - HTMRangeValueVector htmRangeVector = htm->convexHull(points,hullSteps); + HTMRangeValueVector htmRangeVector = htm->convexHull(points,hullSteps,true); // Compress result - cout << dec << "a3000 hrv.size: " << htmRangeVector.size() << endl << flush; + // cout << dec << "a3000 hrv.size: " << htmRangeVector.size() << endl << flush; for(int i=0; i < htmRangeVector.size(); ++i) { uint64 lo = ValueFromHtmID(htmRangeVector[i].lo); // TODO Should this be a function? cover.push_back(lo); uint64 hi; - if( htmRangeVector[i].lo != htmRangeVector[i].lo ) { + if( htmRangeVector[i].lo != htmRangeVector[i].hi ) { hi = sTerminator(ValueFromHtmID(htmRangeVector[i].hi)); cover.push_back(hi); } } - cout << dec << "a4000" << endl << flush; + // cout << dec << "a4000" << endl << flush; delete htm; // TODO Hopefully this will not also delete the index we passed in. return cover; } +STARE_SpatialIntervals STARE::ConvexHull(STARE_ArrayIndexSpatialValues points,int force_resolution_level) { + LatLonDegrees64ValueVector latlon; + for( STARE_ArrayIndexSpatialValues::iterator i=points.begin(); i != points.end(); ++i) { + latlon.push_back(LatLonDegreesFromValue(*i)); + } + return ConvexHull(latlon,force_resolution_level); +} + /* * Return a spatial index object with a given search level. If one does not already exist, construct and memoize. */ @@ -459,6 +556,8 @@ SpatialIndex STARE::getIndex(int resolutionLevel) { /** * Return the legacy htmID value from the spatialStareId. * + * NOTE THIS IGNORES EMBEDDED RESOLUTION + * * Note the htmID precision level needn't have a resolution interpretation, but is more purely geometric. * This is important when calling into the legacy htm foundation and why it's kept private. */ @@ -476,6 +575,7 @@ uint64 STARE::htmIDFromValue(STARE_ArrayIndexSpatialValue spatialStareId, int fo uint64 htmID = rightJustified.getId(); return htmID; } + /** * Return the spatialStareId from the legacy htmID. * @@ -598,26 +698,133 @@ bool terminatorp(STARE_ArrayIndexSpatialValue spatialStareId) { return leftJustifiedWithResolution.terminatorp(); } +STARE_ArrayIndexSpatialValue shiftSpatialIdAtLevel( + STARE_ArrayIndexSpatialValue spatialStareId, + int resolution, + int shiftAmount + ) { +// if( shiftAmount == 0 ) { +// return spatialStareId +// } else { + EmbeddedLevelNameEncoding leftJustifiedWithResolution; + leftJustifiedWithResolution.setIdFromSciDBLeftJustifiedFormat(spatialStareId); + + if( shiftAmount >= 0) { + leftJustifiedWithResolution.setId( + leftJustifiedWithResolution.increment(leftJustifiedWithResolution.getId(), resolution, shiftAmount)); + } else { + leftJustifiedWithResolution.setId( + leftJustifiedWithResolution.decrement(leftJustifiedWithResolution.getId(), resolution, -shiftAmount)); + } + return leftJustifiedWithResolution.getSciDBLeftJustifiedFormat(); +// } +} + +uint64 spatialLevelMask() { + EmbeddedLevelNameEncoding leftJustified; + return leftJustified.levelMaskSciDB; +} + +STARE_ArrayIndexSpatialValues expandInterval(STARE_SpatialIntervals interval, int64 force_resolution) { + // STARE_SpatialIntervals interval should just be one interval, i.e. a value or value+terminator. + DIAGOUT1(cout << endl << dec << 200 << endl << flush;) + STARE_ArrayIndexSpatialValue siv0 = interval[0]; + EmbeddedLevelNameEncoding leftJustified; + DIAGOUT1(cout << dec << 220 << endl << flush;) + uint64 return_resolution = siv0 & leftJustified.levelMaskSciDB; + DIAGOUT1(cout << dec << 225 << " " << setw(16) << setfill('0') << hex << siv0 << dec << endl << flush;) + if( force_resolution > -1 ) { + siv0 = ( siv0 & ~leftJustified.levelMaskSciDB ) | force_resolution; + return_resolution = force_resolution; + } + DIAGOUT1(cout << dec << 229 << " f & resolution: " << dec << force_resolution << " " << return_resolution << endl << flush;) + DIAGOUT1(cout << dec << 230 << " " << setw(16) << setfill('0') << hex << siv0 << dec << endl << flush;) + leftJustified.setIdFromSciDBLeftJustifiedFormat(siv0); + // cout << dec << 235 << endl << flush; + STARE_ArrayIndexSpatialValue siv_term; + if( interval.size() > 1 ) { + siv_term = interval[1]; + } else { + siv_term = leftJustified.getSciDBTerminatorLeftJustifiedFormat(); + } + DIAGOUT1(cout << dec << 239 << " " << setw(16) << setfill('0') << hex << siv_term << dec << endl << flush;) + DIAGOUT1(cout << endl << dec << 240 << endl << flush;) + uint64 one_mask_to_resolution, one_at_resolution; + leftJustified.SciDBincrement_LevelToMaskDelta(siv0 & leftJustified.levelMaskSciDB,one_mask_to_resolution,one_at_resolution); + // cout << dec << 242 << endl << flush; + + uint64 delta = ((siv_term+1)-(siv0 & ~leftJustified.levelMaskSciDB)); + + DIAGOUT1(cout << endl;) + uint64 one = 1; + DIAGOUT1(cout << dec << 243 << " " << setw(16) << setfill('0') << hex << (one << (63-3-2*return_resolution)) << dec << endl << flush;) + DIAGOUT1(cout << dec << 244 << " " << setw(16) << setfill('0') << hex << leftJustified.getSciDBTerminatorLeftJustifiedFormat() << dec << endl << flush;) + DIAGOUT1(cout << dec << 245 << " " << setw(16) << setfill('0') << hex << siv_term << dec << endl << flush;) + DIAGOUT1(cout << dec << 245 << " " << setw(16) << setfill('0') << hex << siv0 << dec << endl << flush;) + DIAGOUT1(cout << dec << 245 << " " << setw(16) << setfill('0') << hex << delta << " " << dec << (delta/one_at_resolution) << endl << flush;) + + // Give as much rope as needed. + DIAGOUT1(cout << dec << 246 << " " << setw(16) << setfill('0') << hex << siv0 << dec << endl << flush;) + siv0 = (siv0 & ~one_mask_to_resolution) | return_resolution; + DIAGOUT1(cout << dec << 247 << " " << setw(16) << setfill('0') << hex << siv0 << dec << endl << endl << flush;) + DIAGOUT1(cout << dec << 247 << " " << setw(16) << setfill('0') << hex << one_at_resolution << dec << endl << endl << flush;) + STARE_ArrayIndexSpatialValues expanded_interval; + while( siv0 < siv_term ) { + DIAGOUT1(cout << dec << 249 << " " << setw(16) << setfill('0') << hex << siv0 << dec << endl << flush;) + expanded_interval.push_back(siv0); + siv0 += one_at_resolution; + } + DIAGOUT1(cout << dec << 250 << endl << endl << flush;) + return expanded_interval; +} + /** - * Construct a range object (spatial region) from a vector of intervals. - * - * TODO: Have a SpatialRange class instead of an HstmRange? + * TODO Fix expandIntervals... */ -HstmRange SpatialRangeFromSpatialIntervals(STARE_SpatialIntervals intervals) { - HstmRange range; +STARE_ArrayIndexSpatialValues expandIntervals(STARE_SpatialIntervals intervals, int64 force_resolution) { + STARE_ArrayIndexSpatialValues expanded_values; EmbeddedLevelNameEncoding leftJustified; - for(auto i0=intervals.begin(); i0 != intervals.end(); ++i0) { - leftJustified.setIdFromSciDBLeftJustifiedFormat(*i0); - uint64 a = leftJustified.getId(), b = a; - auto i1 = (i0+1); - if(terminatorp(*i1)) { - leftJustified.setIdFromSciDBLeftJustifiedFormat(*i1); - b = leftJustified.getId(); - ++i0; // Skip to next + int i=0; + while( i < intervals.size() ) { + // cout << dec << 100 << endl << flush; + STARE_ArrayIndexSpatialValue siv0, siv1; + STARE_SpatialIntervals interval; + siv0 = intervals[i]; + interval.push_back(siv0); + // cout << dec << 110 << endl << flush; + ++i; + if( i < intervals.size() ) { + // peek + // cout << dec << 120 << endl << flush; + siv1 = intervals[i]; + // cout << dec << 120 << " " << i << " " << setw(16) << setfill('0') << hex << siv1 << dec << endl << flush; + if( (siv1 & leftJustified.levelMaskSciDB) == leftJustified.levelMaskSciDB ) { + // cout << dec << 121 << " " << i << endl << flush; + interval.push_back(siv1); + ++i; + } } - range.addRange(a,b); + // cout << dec << 130 << " " << i << endl << flush; + STARE_SpatialIntervals expandOne = expandInterval(interval,force_resolution); + // cout << dec << 140 << " " << i << endl << flush; + for(int j=0; j < expandOne.size(); ++j) { +// cout << dec << 142 << " " << i << " " << j +// << setw(16) << setfill('0') << hex << expandOne[j] << dec +// << endl << flush; + STARE_ArrayIndexSpatialValues_insert( expanded_values, expandOne[j] ); +// cout << dec << 143 << " " << i << " " << j << endl << flush; + } +// cout << dec << 150 << " " << i << endl << flush; } - return range; +// cout << dec << 160 << " " << i << endl << flush; + return expanded_values; +} + +STARE_SpatialIntervals spatialIntervalFromHtmIDKeyPair(KeyPair kp) { + STARE_SpatialIntervals interval; + interval.push_back(EmbeddedLevelNameEncoding(BitShiftNameEncoding(kp.lo).leftJustifiedId()).getSciDBLeftJustifiedFormat()); + interval.push_back(EmbeddedLevelNameEncoding(BitShiftNameEncoding(kp.hi).leftJustifiedId()).getSciDBTerminatorLeftJustifiedFormat()); + return interval; } diff --git a/src/SpatialConstraint.cpp b/src/SpatialConstraint.cpp index 7adcfd1..e9da5fd 100644 --- a/src/SpatialConstraint.cpp +++ b/src/SpatialConstraint.cpp @@ -60,6 +60,14 @@ SpatialConstraint::SpatialConstraint(SpatialVector a, float64 d) : if(d_ >= gEpsilon) sign_ = pOS; } +/////////////CONSTRUCTOR Zero Constraint through two points////////////////////////////////// +// +SpatialConstraint::SpatialConstraint(SpatialVector v0, SpatialVector v1) +{ + a_ = v0 ^ v1; a_.normalize(); d_ = 0; s_ = acos(d_); +} + + /////////////COPY CONSTRUCTOR///////////////////////////// // SpatialConstraint::SpatialConstraint(const SpatialConstraint & old) : diff --git a/src/SpatialIndex.cpp b/src/SpatialIndex.cpp index 8d68780..c2ab035 100644 --- a/src/SpatialIndex.cpp +++ b/src/SpatialIndex.cpp @@ -15,6 +15,7 @@ //# Jul 25, 2002 : Gyorgy Fekete -- Added pointById() //# +#include #include #include @@ -1314,6 +1315,15 @@ SpatialIndex::NeighborsAcrossVerticesFromEdges( workspace[jw++] = q8; } +#ifndef DIAG +#define DIAG1(expr) +#define DIAGOUT2(out,expr) +#define DIAGOUTDELTA(out,a,b) +#else +#define DIAG1(expr) expr; +#define DIAGOUT2(out,expr) out << expr; +#define DIAGOUTDELTA(out,a,b) {SpatialVector delta_ = a-b; cout << delta_.length() << " ";} +#endif //////////////////IDBYPOINT//////////////////////////////////////////////// /** Find a leaf node where a vector points to @@ -1334,24 +1344,158 @@ SpatialIndex::idByPoint(SpatialVector & v) const { uint64 index; uint64 ID; + bool verbose = false; + bool nudge = false; + + /**/ // start with the 8 root triangles, find the one which v points to for(index=1; index <=8; index++) { - if( (V(0) ^ V(1)) * v < -gEpsilon) continue; - if( (V(1) ^ V(2)) * v < -gEpsilon) continue; - if( (V(2) ^ V(0)) * v < -gEpsilon) continue; - break; + if(isInsideBarycentric(v,V(0),V(1),V(2),verbose)) break; } + DIAGOUT2(cout,"0 preamble i: " << index << endl << flush;); + nudge = index == 9; + // loop through matching child until leaves are reached while(ICHILD(0)!=0) { uint64 oldindex = index; for(size_t i = 0; i < 4; i++) { index = nodes_[oldindex].childID_[i]; + if(isInsideBarycentric(v,V(0),V(1),V(2),verbose)) break; + } + } + DIAGOUT2(cout,"1 preamble i: " << index << endl << flush;); + int index_barycentric = index; + /**/ + + // if(index == 9) { + + /* THE OLD WAY + * TODO Why keep the legacy point-finding behavior? Right thing to do would be to examine where barycentric and old way differ. + */ + + SpatialVector dc; + + float64 dcs[8]; + for(index=1; index <=8; index++) { + dc = V(0)+V(1)+V(2); dc.normalize(); dc = dc - v; + dcs[index-1] = dc.length(); + DIAGOUT2(cout,dec << index << " index,dc " << setprecision(16) << dcs[index-1] << endl << flush;); + } + + int index_dcs_sort[8]; + index_dcs_sort[0]=0; + + for(int i = 1; i < 8; ++i) { + int icmp = 0; + while( icmp < i && (dcs[i] >= dcs[index_dcs_sort[icmp]]) ) {++icmp; } + if( icmp <= i) { + for( int j=i+1; j>=icmp; --j ) { + index_dcs_sort[j]=index_dcs_sort[j-1]; + } + } + index_dcs_sort[icmp] = i; + } + DIAG1(for(int i=0;i<8;++i) {DIAGOUT2(cout,i << " i,sort,dcs " << index_dcs_sort[i] << " " << dcs[index_dcs_sort[i]] << endl << flush;);}); + + // start with the 8 root triangles, find the one which v points to + for(index=1; index <=8; index++) { + if( (V(0) ^ V(1)) * v < -gEpsilon) continue; + if( (V(1) ^ V(2)) * v < -gEpsilon) continue; + if( (V(2) ^ V(0)) * v < -gEpsilon) continue; + break; + } + + if(nudge) { + // if(false) { + float64 epsilon = 1.0e-13; // TODO inject this dependency, expose as environment variable + // float64 epsilon = 1.0e-17; // TODO inject this dependency, expose as environment variable + SpatialVector ctr = V(0)+V(1)+V(2); ctr.normalize(); + v = (1-epsilon)*v + epsilon*ctr; v.normalize(); + for(index=1; index <=8; index++) { if( (V(0) ^ V(1)) * v < -gEpsilon) continue; if( (V(1) ^ V(2)) * v < -gEpsilon) continue; if( (V(2) ^ V(0)) * v < -gEpsilon) continue; break; } } + + DIAGOUT2(cout,dec << "index " << index << endl << flush;); + float64 dc_start, dc_end; + float64 dc_improvement = 2.0; + int attempt = 0; + int index_tried = index; + int index_last_found = index; + int itry = 0; + do { + ++attempt; + if(attempt == 9) { + stringstream ss; + float64 lat,lon; + v.getLatLonDegrees(lat, lon); + ss << setprecision(16); + ss << "SpatialIndex::idByPoint(sv): Lost Point Failure 1. No convergence." + << " point = " << v + << " index_barycentric = " << index_barycentric + << " index_last_found = " << index_last_found + << " latlon = " << lat << "," << lon + << " nudge " << ( nudge ? "true" : "false" ) + << " dc_improvement = " << dc_improvement + << " dc_start = " << dc_start + << " dc_end = " << dc_end + << endl << flush; + throw SpatialFailure(ss.str().c_str()); + } + if(attempt>1) { + index = index_dcs_sort[itry]+1; + if( index == index_tried ) { + ++itry; + if(itry == 8) { + stringstream ss; + float64 lat,lon; + v.getLatLonDegrees(lat, lon); + ss << setprecision(16); + ss << "SpatialIndex::idByPoint(sv): Lost Point Failure 2. No convergence." + << " point = " << v + << " index_barycentric = " << index_barycentric + << " index_last_found = " << index_last_found + << " latlon = " << lat << "," << lon + << " nudge " << ( nudge ? "true" : "false" ) + << " dc_improvement = " << dc_improvement + << " dc_start = " << dc_start + << " dc_end = " << dc_end + << endl << flush; + throw SpatialFailure(ss.str().c_str()); + } + index = index_dcs_sort[itry]+1; + } + ++itry; + index_tried = index; + } + dc_start = dcs[index-1]; + DIAG1(cout << attempt << " trying " << index << endl << flush;); + // loop through matching child until leaves are reached + int k = 1; + while(ICHILD(0)!=0) { + uint64 oldindex = index; + for(size_t i = 0; i < 4; i++) { + index = nodes_[oldindex].childID_[i]; + if( (V(0) ^ V(1)) * v < -gEpsilon) continue; + if( (V(1) ^ V(2)) * v < -gEpsilon) continue; + if( (V(2) ^ V(0)) * v < -gEpsilon) continue; + break; + } + dc = V(0)+V(1)+V(2); dc.normalize(); dc = dc - v; + DIAG1(cout << dec << k++ << " k,dc " << setprecision(16) << dc.length() << " i: " << index << endl << flush;); + } + dc = V(0)+V(1)+V(2); dc.normalize(); dc = dc - v; + dc_end = dc.length(); + dc_improvement = dc_end/dc_start; + DIAG1(cout << "dc start, end, ratio: " << dc_start << " " << dc_end << " " << dc_improvement << endl << flush;); + index_last_found = index; + } while ( dc_improvement > 0.125 && dc_end > 0.25 && !nudge ); + // } + /**/ + // what if we haven't gotten to build level? // cout << "idbp: 100 " << endl << flush; // return if we have reached maxlevel @@ -1386,9 +1530,18 @@ SpatialIndex::idByPoint(SpatialVector & v) const { // cout << "idbp: 300 maxlevel_-buildlevel_ = " << level << endl << flush; // TODO make this whole routine less ad-hoc + SpatialVector delta; + DIAGOUT2(cout,endl << flush); if(level>0) { int level0 = buildlevel_; while(level--) { + DIAGOUT2(cout,setprecision(16) << dec;); + // cout << level << " level,v... " << v << " " << v0 << " " << v1 << " " << v2 << endl << flush; + DIAGOUT2(cout,level << " level,d0,d1,d2 ";); + DIAGOUTDELTA(cout,v,v0); + DIAGOUTDELTA(cout,v,v1); + DIAGOUTDELTA(cout,v,v2); + DIAGOUT2(cout,endl << flush;); int subTriangleIndex = subTriangleIndexByPoint(v,v0,v1,v2); // cout << "idbp: 350 level = " << level << ", level0 = " << level0++ << ", subTriangleIndex = " << subTriangleIndex << ", name = " << name; name[len++] = '0'+char(subTriangleIndex); @@ -1400,6 +1553,13 @@ SpatialIndex::idByPoint(SpatialVector & v) const { throw SpatialFailure("SpatialIndex::idByPoint::LevelTooLow len < 2"); } } + + DIAGOUT2(cout,level << " level,d0,d1,d2 ";); + DIAGOUTDELTA(cout,v,v0); + DIAGOUTDELTA(cout,v,v1); + DIAGOUTDELTA(cout,v,v2); + DIAGOUT2(cout,endl << flush << level << " level,v... " << v << " " << v0 << " " << v1 << " " << v2 << endl << flush;); + name[len] = '\0'; // cout << "idbp: 400 name = '" << name << "', len = " << len << endl << flush; @@ -1422,6 +1582,10 @@ SpatialIndex::idByPoint(SpatialVector & v) const { return ID; } +#undef DIAG +#undef DIAG1 +#undef DIAGOUT2 +#undef DIAGOUTDELTA uint64 SpatialIndex::indexAtNodeIndex(uint64 nodeIndex) { return nodes_[nodeIndex].index_; @@ -1457,14 +1621,15 @@ int depthOfId(uint64 htmId) { int levelOfId(uint64 htmId) { int i; uint32 size; + // cout << "loid: " << setw(16) << setfill('0') << hex << htmId << dec << endl << flush; // determine index of first set bit for(i = 0; i < IDSIZE; i+=2) { if ( (htmId << i) & IDHIGHBIT ) break; if ( (htmId << i) & IDHIGHBIT2 ) // invalid id - throw SpatialFailure("SpatialIndex:nameById: invalid ID"); + throw SpatialFailure("SpatialIndex:nameById: invalid ID 1"); } if(htmId == 0) - throw SpatialFailure("SpatialIndex:nameById: invalid ID"); + throw SpatialFailure("SpatialIndex:nameById: invalid ID 2"); size=(IDSIZE-i) >> 1; return size-2; } diff --git a/src/SpatialInterface.cpp b/src/SpatialInterface.cpp index b9f9ccb..689be7d 100644 --- a/src/SpatialInterface.cpp +++ b/src/SpatialInterface.cpp @@ -24,7 +24,7 @@ #define MAX_RANGES 100 // #define DIAG -// #define DIAG_OUT cout +#define DIAG_OUT cout // #define DIAG_OUT cerr #ifdef SpatialSGI @@ -40,6 +40,9 @@ extern long long atoll (const char *str); //============================================================== ///////////CONSTRUCTOR/////////////////////// +htmInterface::htmInterface(const SpatialIndex *index) { + index_ = new SpatialIndex(index->maxlevel_, index->buildlevel_, index->rot_); // TODO delete bait? maxlevel is searchlevel, no? +} htmInterface::htmInterface(size_t searchlevel, size_t buildlevel, SpatialRotation rot) : t_(NULL) { index_ = new SpatialIndex(searchlevel, buildlevel, rot); } @@ -387,24 +390,25 @@ htmInterface::convexHull( } const HTMRangeValueVector & -htmInterface::convexHull( LatLonDegrees64ValueVector latlon, size_t steps ) { +htmInterface::convexHull( LatLonDegrees64ValueVector latlon, size_t steps, bool interiorp ) { + hull_interiorp_ = interiorp; polyCorners_.clear(); - cout << " ch " << 2000 << " latlon-size=" << latlon.size() << flush ; - cout << endl; + // cout << " ch " << 2000 << " latlon-size=" << latlon.size() << flush ; + // cout << endl; if (steps == (uint64) -1) { steps = latlon.size(); } else { steps = min(steps,latlon.size()); } for(size_t i = 0; i < steps; i++) { - cout << " " << i << flush; - cout << " ( " << latlon[i].lat << " " << latlon[i].lon << ")"; + // cout << " " << i << flush; + // cout << " ( " << latlon[i].lat << " " << latlon[i].lon << ")"; float64 *x = xyzFromLatLonDegrees(latlon[i].lat,latlon[i].lon); SpatialVector v(x[0],x[1],x[2]); setPolyCorner(v); - cout << endl << flush; + // cout << endl << flush; } - cout << endl << flush << 2100 << endl << flush; + // cout << endl << flush << 2100 << endl << flush; return doHull(); } @@ -478,7 +482,7 @@ htmInterface::doHull() { // dom.convexes_[0].boundingCircle_.write(cout); // dom.write(cout); //%%%%%%%%%%%%%%%%%%%%%%%%%%% - cout << 3999 << endl << flush; + // cout << 3999 << endl << flush; return domain(dom); } @@ -523,6 +527,25 @@ htmInterface::setPolyCorner(SpatialVector &v) { } else if( (polyCorners_[0].c_ ^ polyCorners_[1].c_)*v < 0 ) { polyCorners_.insert(polyCorners_.end(), polyCorners_[1]); // GYF polyCorners_[1] was missing!!!! polyCorners_[1].c_ = v; + } else { + // Nuts. It's zero. Have to think now. + // Which one is in the middle? + SpatialVector ab = polyCorners_[0].c_ ^ polyCorners_[1].c_; + SpatialVector av = polyCorners_[0].c_ ^ v; + SpatialVector vb = v ^ polyCorners_[1].c_; + float64 dot_ab_av = ab*av; + float64 dot_ab_vb = ab*vb; + + if( dot_ab_av*dot_ab_vb < 0 ) { + // v is not in the middle + if( dot_ab_av > 0 ) { + // b is in the middle + polyCorners_[1].c_ = v; + } else { + // a is in the middle + polyCorners_[0].c_ = v; + } + } // dot_ab_av*dot_ab_vb == 0 is an error and > 0 means v is in the middle. } } else { // @@ -552,6 +575,7 @@ htmInterface::setPolyCorner(SpatialVector &v) { // if( (polyCorners_[i].c_ ^ polyCorners_[i+1==len ? 0 : i+1].c_)*v > tol2 ) { float64 delta = (polyCorners_[i].c_ ^ polyCorners_[i+1==len ? 0 : i+1].c_)*v; +// #define DIAG #ifdef DIAG DIAG_OUT << i << " i," @@ -698,31 +722,34 @@ const HTMRangeValueVector & htmInterface::domain( SpatialDomain & domain ) { HtmRange htmRange; - cout << 4000 << endl << flush; +// cout << 4000 << endl << flush; +// cout << 4001 << " hull_interiorp_ " << hull_interiorp_ << endl << flush; Key gapsize; // SpatialIndex idx(20, 5); // domain.setOlevel(20); - domain.intersect(index_, &htmRange, false); + // domain.intersect(index_, &htmRange, false); + // domain.intersect(index_, &htmRange, true); + domain.intersect(index_, &htmRange, hull_interiorp_); - cout << 4100 << endl << flush; +// cout << 4100 << endl << flush; // gapsize = htmRange.bestgap(MAX_RANGES); // htmRange.defrag(gapsize); - cout << 4200 << endl << flush; -// htmRange.defrag(); +// cout << 4200 << endl << flush; + htmRange.defrag(); // DONT FORGET to: get best gap and defrag htmRange.defrag(bestgap); - cout << 4300 << endl << flush; +// cout << 4300 << endl << flush; // Construct the valuevector... fillValueVec( htmRange, range_); - cout << 4997 << endl << flush; +// cout << 4997 << endl << flush; htmRange.reset(); - cout << 4998 << endl << flush; +// cout << 4998 << endl << flush; htmRange.purge(); - cout << 4999 << endl << flush; +// cout << 4999 << endl << flush; return range_; } diff --git a/src/SpatialRange.cpp b/src/SpatialRange.cpp new file mode 100644 index 0000000..7c0b316 --- /dev/null +++ b/src/SpatialRange.cpp @@ -0,0 +1,137 @@ +/* + * SpatialRange.cpp + * + * Created on: Sep 8, 2019 + * Author: mrilee + * + * Copyright (C) 2019 Rilee Systems Technologies LLC + */ + +#include "SpatialRange.h" + + +SpatialRange::SpatialRange() { + this->range = new HstmRange; +} + +SpatialRange::SpatialRange(STARE_SpatialIntervals intervals) { + this->range = new HstmRange; + this->addSpatialIntervals(intervals); +} + +SpatialRange::~SpatialRange() { + delete this->range; // TODO mlr. Um. Dangerous if range is from elsewhere? Maybe use some sort of smart pointer? +} + +void SpatialRange::addSpatialRange(SpatialRange range) { + this->range->addRange(range.range); +} + +/** + * Construct a range object (spatial region) from a vector of intervals. + * + * TODO: Have a SpatialRange class instead of an HstmRange? + */ +// #define DIAG +void SpatialRange::addSpatialIntervals(STARE_SpatialIntervals intervals) { + EmbeddedLevelNameEncoding leftJustified; + for(auto i0=intervals.begin(); i0 != intervals.end(); ++i0) { + leftJustified.setIdFromSciDBLeftJustifiedFormat(*i0); + uint64 a = leftJustified.getId(), b = a; + auto i1 = (i0+1); + if( i1 <= intervals.end() ) { + if(terminatorp(*i1)) { + leftJustified.setIdFromSciDBLeftJustifiedFormat(*i1); + b = leftJustified.getId(); + ++i0; // Skip to next + } + } + +#ifdef DIAG + cout << "sr::addsi " + << setw(16) << setfill('0') << hex << a << " " + << setw(20) << dec << a << " " + << setw(16) << setfill('0') << hex << b << " " + << setw(20) << dec << b << " " + << endl << flush; +#endif + this->range->addRange(a,b); +#ifdef DIAG + cout << "sr::addsi nr = " << this->range->range->nranges() << endl << flush; +#endif + } +} + +// #define DIAG +int SpatialRange::getNextSpatialInterval(STARE_SpatialIntervals &interval) { + KeyPair kp(-1,-2); + int istat = this->getNext(kp); +#ifdef DIAG + cout << "sr::gnsi istat = " << istat << ", kp = " << kp.lo << ", " << kp.hi << endl << flush; +#endif + if( istat > 0 ) { + EmbeddedLevelNameEncoding leftJustified; + leftJustified.setId(kp.lo); + interval.push_back(leftJustified.getSciDBLeftJustifiedFormat()); + if( kp.lo != kp.hi ) { + STARE_ArrayIndexSpatialValue term_lo = leftJustified.getSciDBTerminatorLeftJustifiedFormat(); + leftJustified.setId(kp.hi); // Note: hi should be a terminator already. + STARE_ArrayIndexSpatialValue term_hi = leftJustified.getSciDBLeftJustifiedFormat(); + +#ifdef DIAG + cout << "lo,term_lo,term_hi: " + << setw(16) << setfill('0') << hex + << kp.lo << "," + << setw(16) << setfill('0') << hex + << term_lo << "," + << setw(16) << setfill('0') << hex + << term_hi << endl << flush << dec; +#endif + + if( term_hi != term_lo ) { + interval.push_back( term_hi ); + } + } + } + return istat; +} +#undef DIAG + +STARE_SpatialIntervals SpatialRange::toSpatialIntervals() { + STARE_SpatialIntervals intervals; + if(this->range) { + this->range->reset(); + while( this->getNextSpatialInterval(intervals) > 0 ); + } + return intervals; +} + +/* + * Odd. The following does not seem to work if we just return the SpatialRange itself. Some of the pointers seem to be either corrupted or eliminated. + */ +SpatialRange *sr_intersect(const SpatialRange&a, const SpatialRange& b, bool compress) { + HstmRange *range = new HstmRange(a.range->range->RangeFromIntersection(b.range->range,compress)); // NOTE mlr Probably about the safest way to inst. SpatialRange. +// #define DIAG +#ifdef DIAG + KeyPair kp; range->reset(); range->getNext(kp); + cout << "sr_i range,r->r,nr " << range << " " << range->range << " " << range->range->nranges() << " : " + << setw(16) << setfill('0') << hex << kp.lo << " " + << setw(16) << setfill('0') << hex << kp.hi << " " + << dec + << endl << flush; + EmbeddedLevelNameEncoding leftJustified; + leftJustified.setId(kp.lo); cout << "kp.lo lj " << setw(16) << setfill('0') << hex << leftJustified.getSciDBLeftJustifiedFormat() << endl << flush; + leftJustified.setId(kp.hi); cout << "kp.hi lj " << setw(16) << setfill('0') << hex << leftJustified.getSciDBLeftJustifiedFormat() << endl << flush; + cout << " r-r-my_los " << hex << range->range->my_los << endl << flush; + cout << dec; +#endif + SpatialRange *sr = new SpatialRange(range); +#ifdef DIAG + cout << "sr-r " << hex << sr->range << endl << flush; + cout << "sr-r-r " << hex << sr->range->range << endl << flush; + cout << "sr-r-r-my_los " << hex << sr->range->range->my_los << endl << flush; +#endif + return sr; +} + + diff --git a/tests/CUTE/CMakeLists.txt b/tests/CUTE/CMakeLists.txt index 4d4b2e5..f8b4089 100644 --- a/tests/CUTE/CMakeLists.txt +++ b/tests/CUTE/CMakeLists.txt @@ -1,15 +1,17 @@ set ( TestSrcFiles - Test.cpp + Test.h + EmbeddedLevelNameEncoding_test.cpp SpatialFailure_test.cpp SpatialIndex_test.cpp SpatialInterface_test.cpp + SpatialRange_test.cpp SpatialRotation_test.cpp SpatialVector_test.cpp STARE_test.cpp TemporalIndex_test.cpp - Test.h + Test.cpp ) # SET( CMAKE_EXE_LINKER_FLAGS "${CMAKE_EXE_LINKER_FLAGS} -I${PROJECT_SOURCE_DIR}/include" ) diff --git a/tests/CUTE/EmbeddedLevelNameEncoding_test.cpp b/tests/CUTE/EmbeddedLevelNameEncoding_test.cpp new file mode 100644 index 0000000..f4f2e87 --- /dev/null +++ b/tests/CUTE/EmbeddedLevelNameEncoding_test.cpp @@ -0,0 +1,100 @@ +/* + * EmbeddedLevelNameEncoding_test.cpp + * + * Created on: Sep 9, 2019 + * Author: mrilee + * + * Copyright (C) 2019 Rilee Systems Technologies LLC + */ + +#include "Test.h" + +void EmbeddedLevelNameEncoding_test() { + STARE index; + STARE_ArrayIndexSpatialValue sid = index.ValueFromLatLonDegrees( 30, 30, 8 ); + + EmbeddedLevelNameEncoding leftJustified; + + leftJustified.setIdFromSciDBLeftJustifiedFormat(sid); +#ifdef DIAG + cout << "sid tp? " << leftJustified.terminatorp() << endl << flush; +#endif + + STARE_ArrayIndexSpatialValue term_sid1 = leftJustified.getSciDBTerminatorLeftJustifiedFormat(); + leftJustified.setIdFromSciDBLeftJustifiedFormat(term_sid1); +#ifdef DIAG + cout << "term_sid1 tp? " << leftJustified.terminatorp() << endl << flush; +#endif + + STARE_ArrayIndexSpatialValue term_sid2 = leftJustified.getSciDBTerminatorLeftJustifiedFormat(); + +#ifdef DIAG + cout << "test: id,mask,level " + << setw(16) << setfill('0') << hex + << leftJustified.getId() + << " " << dec << leftJustified.levelMask + << " " << (leftJustified.getId() & leftJustified.levelMask) + << endl << flush; + cout << "testS: id,mask,level " + << setw(16) << setfill('0') << hex + << leftJustified.getId() + << " " << dec << leftJustified.levelMaskSciDB + << " " << (leftJustified.getId() & leftJustified.levelMaskSciDB) + << endl << flush; + cout << endl << flush; + +#define IDOUT(p,m,s) p << m << " " << setw(16) << setfill(' ') << hex << s << dec << endl << flush; + IDOUT(cout,"sid ", sid); + IDOUT(cout,"term_sid1 ", term_sid1); + IDOUT(cout,"term_sid1*", leftJustified.getId()); + IDOUT(cout,"term_sid1*", leftJustified.getSciDBLeftJustifiedFormat()); + IDOUT(cout,"term_sid1t", leftJustified.getSciDBTerminatorLeftJustifiedFormat()); + IDOUT(cout,"term_sid2 ", term_sid2); +#undef IDOUT +#endif + + ASSERT(leftJustified.terminatorp(term_sid1)); + ASSERT_EQUALM("getSciDBTerminatorLeftJustifiedFormat idempotent?",term_sid1,term_sid2); + + if(true) { + STARE_ArrayIndexSpatialValue siv0 = 0x2324000000000005; + EmbeddedLevelNameEncoding lj; lj.setIdFromSciDBLeftJustifiedFormat(siv0); + STARE_ArrayIndexSpatialValue siv0_term = lj.getSciDBTerminatorLeftJustifiedFormat(); + + STARE_ArrayIndexSpatialValue not_siv1 = lj.increment(lj.getId(), lj.getLevel(), 1); + EmbeddedLevelNameEncoding lj1(not_siv1); + STARE_ArrayIndexSpatialValue siv1 = lj1.getSciDBLeftJustifiedFormat(); + STARE_ArrayIndexSpatialValue sivt_term = 0x2327ffffffffffff; + + int level = lj.getLevel(); + uint64 one_mask_to_level, one_at_level; + lj.SciDBincrement_LevelToMaskDelta(level, one_mask_to_level, one_at_level); + STARE_ArrayIndexSpatialValue siv2 = siv0 + one_at_level; + STARE_ArrayIndexSpatialValue siv2_term = siv2 | one_mask_to_level; + + STARE_ArrayIndexSpatialValue siv0_tchk = (siv0 | one_mask_to_level); + +#define IDOUT(p,m,s) +// #define IDOUT(p,m,s) p << m << " " << setw(16) << setfill(' ') << hex << s << dec << endl << flush; + IDOUT(cout,"siv0 ",siv0); + IDOUT(cout,"siv0_term ",siv0_term); + IDOUT(cout,"siv0_tchk ",siv0_tchk); + IDOUT(cout,"not_siv1 ",not_siv1); + IDOUT(cout,"siv1 ",siv1); + IDOUT(cout,"sivt_term ",sivt_term); + IDOUT(cout,"siv2 ",siv2); + IDOUT(cout,"siv2_term ",siv2_term); + // IDOUT(cout,"lj id ",lj.getId()); +#undef IDOUT + ASSERT_EQUAL( siv0 , 0x2324000000000005 ); + ASSERT_EQUAL( siv0_term , 0x2325ffffffffffff ); + ASSERT_EQUAL( siv0_tchk , 0x2325ffffffffffff ); + ASSERT_EQUAL( not_siv1 , 0xc64c000000000005 ); + ASSERT_EQUAL( siv1 , 0x2326000000000005 ); + ASSERT_EQUAL( sivt_term , 0x2327ffffffffffff ); + ASSERT_EQUAL( siv2 , 0x2326000000000005 ); + ASSERT_EQUAL( siv2_term , 0x2327ffffffffffff ); + } + + +} diff --git a/tests/CUTE/STARE_test.cpp b/tests/CUTE/STARE_test.cpp index 8522420..3b93c2c 100644 --- a/tests/CUTE/STARE_test.cpp +++ b/tests/CUTE/STARE_test.cpp @@ -12,6 +12,16 @@ #include "SpatialInterface.h" +#ifndef DIAG +#define DIAG1(expr) +#define DIAGOUT2(out,expr) +#define DIAGOUTDELTA(out,a,b) +#else +#define DIAG1(expr) expr; +#define DIAGOUT2(out,expr) out << expr; +#define DIAGOUTDELTA(out,a,b) {SpatialVector delta_ = a-b; cout << delta_.length() << " ";} +#endif + void STARE_test() { SpatialVector @@ -908,6 +918,391 @@ void STARE_test() { // cout << " idx=0, latlon0: " << latlon0.lat << " " << latlon0.lon << endl << flush; // } -// FAIL(); + + if(false) { + LatLonDegrees64 latlon(32.5735,-100.05); + STARE index2; + int resolution = 1; + STARE_ArrayIndexSpatialValue idx = index2.ValueFromLatLonDegrees(latlon.lat,latlon.lon,resolution); + + idx = shiftSpatialIdAtLevel( idx, resolution, 1 ); + + cout << endl << "base..." << endl << flush; + cout + << setprecision(18) + << setw(23) + << " idx = 0x" << hex << idx << dec + << scientific + << " latlon = " + << latlon.lat << "," << latlon.lon + << dec << " rLevel = " << resolution + << endl << flush; + + + cout << endl << "increasing resolution" << endl << flush; + for(int resolution_ = 0; resolution_ < 12; ++ resolution_) { + STARE_ArrayIndexSpatialValue idx1 = shiftSpatialIdAtLevel( idx, resolution_, 1 ); + LatLonDegrees64 latlon1 = index2.LatLonDegreesFromValue(idx1); + int resolution1 = idx1 & spatialLevelMask(); + + cout + << setprecision(18) + << setw(23) + << " idx = 0x" << hex << idx1 << dec + << scientific + << " latlon = " + << latlon1.lat << "," << latlon1.lon + << dec << " rLevel = " << resolution1 + << " resLevel = " << resolution_ + << endl << flush; + } + + cout << endl << "decreasing resolution" << endl << flush; + for(int resolution_ = 0; resolution_ < 12; ++ resolution_) { + STARE_ArrayIndexSpatialValue idx1 = shiftSpatialIdAtLevel( idx, resolution_, -1 ); + LatLonDegrees64 latlon1 = index2.LatLonDegreesFromValue(idx1); + int resolution1 = idx1 & spatialLevelMask(); + + cout + << setprecision(18) + << setw(23) + << " idx = 0x" << hex << idx1 << dec + << scientific + << " latlon = " + << latlon1.lat << "," << latlon1.lon + << dec << " rLevel = " << resolution1 + << " resLevel = " << resolution_ + << endl << flush; + } + + { + cout << endl << "incrementing at resolution" << endl << flush; + uint64 resolution_ = 0; + STARE_ArrayIndexSpatialValue idx1 = 0x0000000000000001; + stringstream ss; + + for( resolution_ = 0; resolution_ < 8; ++resolution_) { + for(uint64 itmp = 0; itmp < 7; ++itmp) { + + // cout << resolution_ << "," << itmp << " res,i" << endl << flush; + + try { + idx1 = shiftSpatialIdAtLevel( 0x0000000000000000+resolution_, resolution_, itmp ); + // idx1 = shiftSpatialIdAtLevel( 0ul+resolution_, resolution_, itmp ); + } catch (const SpatialException & e ) { + cout << "Exception: " << dec << e.what() << endl << flush; + FAIL(); + } + + LatLonDegrees64 latlon1 = index2.LatLonDegreesFromValue(idx1); + + int resolution1 = idx1 & spatialLevelMask(); + + ss.clear(); ss.str(string()); + ss << "Increment by " << itmp << " at resolution level " << resolution_; + + cout + << " idx = 0x" << setw(16) << hex << idx1 << dec + << " cmp = 0x" << setw(16) << hex << (itmp << 59) << dec + << scientific << setprecision(18) << setw(23) + << " latlon = " + << setfill(' ') + << latlon1.lat << "," << latlon1.lon + << dec << " rLevel = " << resolution1 + << " resLevel = " << resolution_ + << endl << flush; + + ASSERT_EQUALM(ss.str().c_str(),((itmp << (59 - 2*resolution_)) || resolution_),idx1); + + } + } + } + } + +// #define DIAG +#ifndef DIAG +#define DIAGOUT2(p,m) +#define SIVOUT(m,siv) +#define SIVSOUT(p,m,v) +#else +#define DIAGOUT2(p,m) p << m; +#define SIVOUT(m,siv) cout << m << " " << setw(16) << setfill('0') << hex << siv << dec << endl << flush; +#define SIVSOUT(p,m,v) { p << m << " "; for(int l=0; lpolyCorners_[3].c_,v,1.0e-12); delete htm; } + + if( true ) { + + // cout << endl << "----" << endl << endl << endl << flush; + + // Try a grid... + int nLat = 11,nLon = 11, nLatLon=nLat*nLon; + float64 lat[nLat], lon[nLon]; + SpatialVector vs[nLatLon], cs[nLatLon], vs1[nLatLon]; + for( int i = 0; i < nLon; ++i ) { + lon[i] = i*5.0; + } + for( int j = 0; j < nLat; ++j ) { + lat[j] = j*5.0; + } + int k = 0; + for( int i = 0; i < nLon; ++i ) { + for( int j = 0; j < nLat; ++j ) { + vs[k++].setLatLonDegrees(lat[j], lon[i]); + } + } + + STARE index; + htmInterface *htm; + int resolution_level = 4; + htm = new htmInterface( + index.getIndex(resolution_level).getMaxlevel(), + index.getIndex(resolution_level).getBuildLevel(), + index.getIndex(resolution_level).getRotation()); + htm->polyCorners_.clear(); + SpatialVector v; + + k = 0; + for( int i=0; isetPolyCorner(vs[k]); + // cout << k << " k, hpc size " << htm->polyCorners_.size() << endl << flush; + float64 lat_,lon_; vs[k].getLatLonDegrees(lat_,lon_); + // cout << k << " k, try " + // << lat_ << " " << lon_ << endl << flush; + for(int l=0; lpolyCorners_.size(); ++l) { + htm->polyCorners_[l].c_.getLatLonDegrees(lat_,lon_); + // cout << l << " l, hpc " + // << lat_ << " " << lon_ << endl << flush; + } + // cout << "-----" << endl << flush; + ++k; + // if(k>5) exit(1); + } + } + + delete htm; + + } } } diff --git a/tests/CUTE/SpatialRange_test.cpp b/tests/CUTE/SpatialRange_test.cpp new file mode 100644 index 0000000..444acec --- /dev/null +++ b/tests/CUTE/SpatialRange_test.cpp @@ -0,0 +1,282 @@ +/* + * SpatialRange_test.cpp + * + * Created on: Sep 8, 2019 + * Author: mrilee + * + * Copyright (C) 2019 Rilee Systems Technologies LLC + */ + +#include "Test.h" + +// #define DIAG + +void SpatialRange_test () { + + STARE index; + STARE_ArrayIndexSpatialValue siv = index.ValueFromLatLonDegrees( 30, 30, 8 ); + STARE_SpatialIntervals sis; sis.push_back(siv); + + SpatialRange range(sis); + STARE_SpatialIntervals sis_out; + + sis_out = range.toSpatialIntervals(); + +// #define DIAG +#ifdef DIAG + for(int i=0; i < sis.size(); ++i ) { + cout << i << " i,si: " << setw(16) << setfill('0') << hex << sis[i] << dec << endl << flush; + } + for(int i=0; i < sis_out.size(); ++i ) { + cout << i << " i,so: " << setw(16) << setfill('0') << hex << sis_out[i] << dec << endl << flush; + } +#endif + + ASSERT_EQUAL(1,sis_out.size()); + ASSERT_EQUAL(sis[0],sis_out[0]); + + // sids.push_back(index.ValueFromLatLonDegrees( 45, 45, 8)); + + siv = index.ValueFromLatLonDegrees( 15, 15, 8 ); + + EmbeddedLevelNameEncoding leftJustified; + leftJustified.setIdFromSciDBLeftJustifiedFormat(siv); + sis.push_back(leftJustified.getSciDBTerminatorLeftJustifiedFormat()); + +#ifdef DIAG + for(int i=0; i < sis.size(); ++i ) { + cout << i << " i,si: " << setw(16) << setfill('0') << hex << sis[i] << dec << endl << flush; + } +// for(int i=0; i < sis_out.size(); ++i ) { +// cout << i << " i,so: " << setw(16) << setfill('0') << hex << sis_out[i] << dec << endl << flush; +// } +#endif + + // cout << 50 << endl << flush; + + range.purge(); + range.addSpatialIntervals(sis); +// cout << 100 +// << " nr: " << range.range->range->nranges() +// << endl << flush; + sis_out = range.toSpatialIntervals(); + +#ifdef DIAG + cout << 110 << " sis_out.s: " << sis_out.size() << endl << flush; + for(int i=0; irange) << dec << endl << flush; + cout << 1991 << " dR.r->r " << hex << (deltaRange->range->range) << dec << endl << flush; + cout << 1991 << " dR.r->r->my_los " << hex << (deltaRange->range->range->my_los) << dec << endl << flush; +#endif + if(!deltaRange->range->range) { cout << "Error deltaRange is null!" << endl << flush; } + // cout << 200 << " dR nR = " << deltaRange->range->range->nranges() << endl << flush; + STARE_SpatialIntervals deltaSis = deltaRange->toSpatialIntervals(); + // cout << 300 << endl << flush; + SISOUT("deltaSis",deltaSis); + // cout << 400 << endl << flush; + + ASSERT_EQUAL(sis1[0],deltaSis[0]); + ASSERT_EQUAL(sis0[1],deltaSis[1]); + +#undef SIVOUT +#undef SISOUT + + } + + + if(true) { + // TODO Fix the following to use the new 0/1 variables defined above. +// #define DIAG +#ifndef DIAG +#define SIVOUT(m,siv) +#define SISOUT(m,sis) +#else +#define SIVOUT(m,siv) cout << m << " " << setw(16) << setfill('0') << hex << siv << dec << endl << flush; +#define SISOUT(m,sis) { cout << endl << m << endl; for(int i=0; i < sis.size(); ++i ) { cout << i << " i,si: " << setw(16) << setfill('0') << hex << sis[i] << dec << endl << flush; }} +#endif + // #define DIAGOUT + // #define DIAGOUT(m) {cout << m << endl << flush;} + // Note level is 8. + uint64 one_mask_to_level, one_at_level; + leftJustified.SciDBincrement_LevelToMaskDelta(8,one_mask_to_level, one_at_level); + + STARE_SpatialIntervals sis0; + SpatialRange range0; + STARE_ArrayIndexSpatialValue siv0 = 0x00000000000000008; + sis0.push_back(siv0); + siv0+=6*one_at_level; + leftJustified.setIdFromSciDBLeftJustifiedFormat(siv0); // TODO Rename this STARE, or provide a STARE interface for this. + sis0.push_back(leftJustified.getSciDBTerminatorLeftJustifiedFormat()); + range0.addSpatialIntervals(sis0); + SISOUT("sis0",sis0); + + STARE_SpatialIntervals sis1; + SpatialRange range1; + STARE_ArrayIndexSpatialValue siv1 = 0x00000000000000008; + siv1+=3*one_at_level; + siv1 = ( siv1 & ~leftJustified.levelMaskSciDB ) | 10; + sis1.push_back(siv1); + siv1+=6*one_at_level; + leftJustified.setIdFromSciDBLeftJustifiedFormat(siv1); + sis1.push_back(leftJustified.getSciDBTerminatorLeftJustifiedFormat()); + range1.addSpatialIntervals(sis1); + SISOUT("sis1",sis1); + + // cout << 100 << endl << flush; + SpatialRange *deltaRange; + try { + // cout << 110 << endl << flush; + deltaRange = range0 & range1; + // deltaRange = sr_intersect(range0, range1); + // cout << 120 << endl << flush; + } catch ( SpatialException e ) { + // cout << 130 << endl << flush; + cout << e.what() << endl << flush; + // cout << 140 << endl << flush; + } + +#ifdef DIAG + cout << 199 << endl << flush; + cout << 1991 << " dR.r " << hex << (deltaRange->range) << dec << endl << flush; + cout << 1991 << " dR.r->r " << hex << (deltaRange->range->range) << dec << endl << flush; + cout << 1991 << " dR.r->r->my_los " << hex << (deltaRange->range->range->my_los) << dec << endl << flush; +#endif + if(!deltaRange->range->range) { cout << "Error deltaRange is null!" << endl << flush; } + // cout << 200 << " dR nR = " << deltaRange->range->range->nranges() << endl << flush; + STARE_SpatialIntervals deltaSis = deltaRange->toSpatialIntervals(); + // cout << 300 << endl << flush; + SISOUT("deltaSis",deltaSis); + // cout << 400 << endl << flush; + + ASSERT_EQUAL(2,deltaSis.size()); + ASSERT_EQUAL(sis1[0],deltaSis[0]); + ASSERT_EQUAL(sis0[1],deltaSis[1]); + } + + // TODO Write many more tests & consider edge cases. + if(true) { + STARE_ArrayIndexSpatialValue siv1[2] = { 0x0000000000000008, 0x000067ffffffffff }; + STARE_ArrayIndexSpatialValue siv2[2] = { 0x000030000000000a, 0x0000907fffffffff }; + STARE_SpatialIntervals sis1(siv1,siv1+2); + STARE_SpatialIntervals sis2(siv2,siv2+2); + SpatialRange r1(sis1), r2(sis2); + SpatialRange *ri = r1 & r2; + // SpatialRange *ri = sr_intersect(r1, r2); + + STARE_SpatialIntervals result = ri->toSpatialIntervals(); + ASSERT_EQUAL(2,result.size()); + ASSERT_EQUAL(0x000030000000000a,result[0]); + ASSERT_EQUAL(0x000067ffffffffff,result[1]); + + ri->range->range->CompressionPass(); + result = ri->toSpatialIntervals(); + ASSERT_EQUAL(4,result.size()); + ASSERT_EQUAL(0x0000300000000008,result[0]); + ASSERT_EQUAL(0x00003fffffffffff,result[1]); + ASSERT_EQUAL(0x0000400000000007,result[2]); + ASSERT_EQUAL(0x0000600000000008,result[3]); + + delete ri; + ri = sr_intersect(r1, r2, true); // Run a compression pass on the range. + result.clear(); + result = ri->toSpatialIntervals(); + ASSERT_EQUAL(4,result.size()); + ASSERT_EQUAL(0x0000300000000008,result[0]); + ASSERT_EQUAL(0x00003fffffffffff,result[1]); + ASSERT_EQUAL(0x0000400000000007,result[2]); + ASSERT_EQUAL(0x0000600000000008,result[3]); + } + + if(false) { + // From PySTARE.cpp + // void _intersect(int64_t* indices1, int len1, int64_t* indices2, int len2, int64_t* intersection, int leni) { + int len1 = 3, len2 = 3, leni = 3; + int64_t + indices1[len1] = {0,0,0}, + // indices2[len2] = {0,0,0}, + intersection[leni]; + + int64_t *indices2 = indices1; + + STARE_SpatialIntervals si1(indices1, indices1+len1), si2(indices2, indices2+len2); + // intersection[0] = 69; + SpatialRange r1(si1), r2(si2); + SpatialRange *ri = r1 & r2; + STARE_SpatialIntervals result = ri->toSpatialIntervals(); + leni = result.size(); + for(int i=0; i