… | |
… | |
17 | memcpy (n, k, ndims * sizeof (double)); |
17 | memcpy (n, k, ndims * sizeof (double)); |
18 | return n; |
18 | return n; |
19 | } |
19 | } |
20 | */ |
20 | */ |
21 | |
21 | |
|
|
22 | double *Object::data() const |
|
|
23 | { |
|
|
24 | double *d = new double [NDIMS]; |
|
|
25 | |
|
|
26 | for (int i = NDIMS; i--; ) |
|
|
27 | d[i] = int2double(k[i]); |
|
|
28 | |
|
|
29 | return d; |
|
|
30 | } |
|
|
31 | |
22 | double Object::distance(const Object& other) const |
32 | double Object::distance(const Object& other) const |
23 | { |
33 | { |
24 | int i; |
34 | if (ACC->distfast) |
25 | double dist = 0; |
|
|
26 | |
|
|
27 | for (i = NDIMS; i--; ) |
|
|
28 | { |
35 | { |
|
|
36 | velem dist = 0; |
|
|
37 | |
|
|
38 | for (int i = NDIMS; i--; ) |
|
|
39 | { |
29 | double d = other.k[i] - k[i]; |
40 | long d = other.k[i] - k[i]; |
30 | |
41 | |
31 | dist += d*d; |
42 | dist += (velem)(d*d); |
|
|
43 | } |
|
|
44 | |
|
|
45 | return dist * ACC->distmul; |
32 | } |
46 | } |
|
|
47 | else |
|
|
48 | { |
|
|
49 | double dist = 0.; |
33 | |
50 | |
|
|
51 | for (int i = NDIMS; i--; ) |
|
|
52 | { |
|
|
53 | double d = int2double(other.k[i]) - int2double(k[i]); |
|
|
54 | |
|
|
55 | dist += (velem)(d*d); |
|
|
56 | } |
|
|
57 | |
34 | return dist; |
58 | return dist; |
|
|
59 | } |
|
|
60 | } |
|
|
61 | |
|
|
62 | |
|
|
63 | Object::Object(double *pkey) |
|
|
64 | { |
|
|
65 | k = new velem[NDIMS]; |
|
|
66 | |
|
|
67 | // discretize the vector |
|
|
68 | for (int i = NDIMS; i--; ) |
|
|
69 | k[i] = double2int (pkey[i]); |
|
|
70 | |
|
|
71 | free (pkey); |
35 | } |
72 | } |
36 | |
73 | |
37 | Object::Object(char *key) |
74 | Object::Object(char *key) |
38 | { |
75 | { |
39 | unsigned char *c = (unsigned char *)key; |
76 | unsigned char *c = (unsigned char *)key; |
40 | |
77 | |
41 | k = new double [NDIMS]; |
78 | k = new velem [NDIMS]; |
42 | |
79 | |
43 | for (int i = 0; i < NDIMS; i++) |
80 | for (int i = 0; i < NDIMS; i++) |
44 | { |
81 | { |
45 | unsigned long elem = 0; |
82 | velem elem = 0; |
46 | |
83 | |
47 | switch (ACC->elemsize) |
84 | switch (ACC->elemsize) |
48 | { |
85 | { |
49 | case 4: elem |= *c++ << 24; |
86 | case 4: elem |= *c++ << 24; |
50 | case 3: elem |= *c++ << 16; |
87 | case 3: elem |= *c++ << 16; |
… | |
… | |
53 | break; |
90 | break; |
54 | default: |
91 | default: |
55 | abort (); |
92 | abort (); |
56 | } |
93 | } |
57 | |
94 | |
58 | k[i] = ((double)elem) * ACC->max / ACC->steps + ACC->min; |
95 | k[i] = elem; |
59 | } |
96 | } |
60 | } |
97 | } |
61 | |
98 | |
62 | void Object::Compress(char *key) |
99 | void Object::Compress(char *key) |
63 | { |
100 | { |
64 | unsigned char *c = (unsigned char *)key; |
101 | unsigned char *c = (unsigned char *)key; |
65 | |
102 | |
66 | for (int i = 0; i < NDIMS; i++) |
103 | for (int i = 0; i < NDIMS; i++) |
67 | { |
104 | { |
68 | unsigned long elem = (unsigned long)floor ((k[i] - ACC->min) * ACC->steps / ACC->max); |
105 | velem elem = k[i]; |
69 | |
106 | |
70 | switch (ACC->elemsize) |
107 | switch (ACC->elemsize) |
71 | { |
108 | { |
72 | case 4: *c++ = (elem >> 24) & 0xff; |
109 | case 4: *c++ = (elem >> 24) & 0xff; |
73 | case 3: *c++ = (elem >> 16) & 0xff; |
110 | case 3: *c++ = (elem >> 16) & 0xff; |
… | |
… | |
80 | } |
117 | } |
81 | } |
118 | } |
82 | |
119 | |
83 | #define SETCUR current_pmt = this |
120 | #define SETCUR current_pmt = this |
84 | |
121 | |
85 | PMT::PMT(int ndims, double min, double max, double steps) |
122 | PMT::PMT(int ndims, double min, double max, double steps, unsigned int pagesize) |
86 | { |
123 | { |
|
|
124 | steps--; |
|
|
125 | |
|
|
126 | assert (min <= 0); |
|
|
127 | |
87 | this->ndims = ndims; |
128 | this->ndims = ndims; |
88 | this->min = min; |
129 | this->min = min; |
89 | this->max = max; |
130 | this->max = max; |
90 | this->steps = steps; |
131 | this->steps = steps; |
91 | |
132 | |
|
|
133 | this->vzero = (velem)floor (- min * steps / max); |
|
|
134 | |
92 | if (steps <= (1<<8)) |
135 | if (steps < (1<<8)) |
93 | elemsize = 1; |
136 | elemsize = 1; |
94 | else if (steps <= (1<<16)) |
137 | else if (steps < (1<<16)) |
95 | elemsize = 2; |
138 | elemsize = 2; |
96 | else if (steps <= (1<<24)) |
139 | else if (steps < (1<<24)) |
97 | elemsize = 3; |
140 | elemsize = 3; |
98 | else |
141 | else |
99 | elemsize = 4; |
142 | elemsize = 4; |
100 | |
143 | |
101 | maxDist = (max - min) * (max - min) * ndims; |
144 | maxDist = (max - min) * (max - min) * ndims; |
102 | |
145 | |
|
|
146 | if (elemsize <= 2 |
|
|
147 | && floor (steps * steps * 0.5 + 1) < velem(1<<31) / velem(ndims)) |
|
|
148 | { |
|
|
149 | distfast = 1; |
|
|
150 | distmul = maxDist / double(steps * steps * ndims); |
|
|
151 | } |
|
|
152 | else |
|
|
153 | distfast = 0; |
|
|
154 | |
103 | SETCUR; |
155 | SETCUR; |
104 | mt = new MT; |
156 | mt = new MT (pagesize); |
105 | } |
157 | } |
106 | |
158 | |
107 | PMT::~PMT() |
159 | PMT::~PMT() |
108 | { |
160 | { |
109 | SETCUR; |
161 | SETCUR; |
|
|
162 | mt->Sync(); |
110 | delete mt; |
163 | delete mt; |
111 | } |
164 | } |
112 | |
165 | |
113 | void PMT::create(const char *path) |
166 | void PMT::create(const char *path) |
114 | { |
167 | { |
… | |
… | |
129 | const MTkey key(o, 0, 0); |
182 | const MTkey key(o, 0, 0); |
130 | const MTentry entry(key, data); |
183 | const MTentry entry(key, data); |
131 | mt->Insert(entry); |
184 | mt->Insert(entry); |
132 | } |
185 | } |
133 | |
186 | |
|
|
187 | void PMT::sync() |
|
|
188 | { |
|
|
189 | SETCUR; |
|
|
190 | mt->Sync(); |
|
|
191 | } |
|
|
192 | |
134 | double PMT::distance(double *k1, double *k2) const |
193 | double PMT::distance(double *k1, double *k2) const |
135 | { |
194 | { |
136 | SETCUR; |
195 | SETCUR; |
137 | Object o1(k1), o2(k2); |
196 | Object o1(k1), o2(k2); |
138 | return o1.distance(o2); |
197 | return o1.distance(o2); |
… | |
… | |
145 | Object o(k); |
204 | Object o(k); |
146 | Pred p(o); |
205 | Pred p(o); |
147 | SimpleQuery q(&p, r); |
206 | SimpleQuery q(&p, r); |
148 | GiSTlist<MTentry *> res = mt->RangeSearch(q); |
207 | GiSTlist<MTentry *> res = mt->RangeSearch(q); |
149 | |
208 | |
150 | while(!res.IsEmpty()) { |
209 | while(!res.IsEmpty()) |
|
|
210 | { |
151 | MTentry *e = res.RemoveFront (); |
211 | MTentry *e = res.RemoveFront (); |
152 | |
212 | double *data = e->Key()->obj.data (); |
153 | add_result(e->Key()->distance, e->Ptr(), e->Key()->obj.data(), ndims); |
213 | add_result(e->Ptr(), data, ndims); |
154 | |
214 | delete data; |
155 | delete e; |
215 | delete e; |
156 | } |
216 | } |
157 | } |
217 | } |
158 | |
218 | |
159 | void PMT::top(double *k, int n) const |
219 | void PMT::top(double *k, int n) const |
160 | { |
220 | { |
161 | abort (); |
221 | SETCUR; |
|
|
222 | Object o(k); |
|
|
223 | Pred p(o); |
|
|
224 | TopQuery q(&p, n); |
|
|
225 | MTentry **res = mt->TopSearch(q); |
|
|
226 | |
|
|
227 | for (int i=0; i < n; i++) |
|
|
228 | { |
|
|
229 | MTentry *e = res[i]; |
|
|
230 | add_result(e->Ptr(), e->Key()->obj.data(), ndims); |
|
|
231 | delete e; |
|
|
232 | } |
162 | } |
233 | } |
|
|
234 | |
|
|
235 | int PMT::maxlevel() const |
|
|
236 | { |
|
|
237 | SETCUR; |
|
|
238 | return mt->MaxLevel(); |
|
|
239 | } |