Skip to content

Commit 4b4fcff

Browse files
bhallerbhaller
committed
switch to binary search for findInterval() for performance
1 parent 0c5c98f commit 4b4fcff

2 files changed

Lines changed: 96 additions & 0 deletions

File tree

‎VERSIONS‎

Lines changed: 1 addition & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -175,6 +175,7 @@ development head (in the master branch):
175175
add xy, xz, yz, and xyz properties to Individual that allow joint access to the x/y/z properties, for easier spatiality operations in complex spatial models (and to support some unit testing additions)
176176
fixed a variety of memory leaks and extended the unit testing for leaks considerably; should be pretty clean in its memory usage now
177177
switch chapters 15 and 16 (nonWF and spatial models) and reorganize/renumber to rationalize the manual organization; many recipe numbers changed
178+
speed up findInterval() by using binary search instead of linear search - O(log n) instead of O(n) to find the interval for one element, with n being the number of intervals (i.e., the length of vec)
178179

179180

180181
version 4.0.1 (Eidos version 3.0.1):

‎eidos/eidos_functions_distributions.cpp‎

Lines changed: 95 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -56,6 +56,9 @@
5656
// (integer)findInterval(numeric x, numeric vec, [logical$ rightmostClosed = F], [logical$ allInside = F])
5757
EidosValue_SP Eidos_ExecuteFunction_findInterval(const std::vector<EidosValue_SP> &p_arguments, __attribute__((unused)) EidosInterpreter &p_interpreter)
5858
{
59+
// Switching to using binary search for this algorithm, but I want to keep the old code around in case it is needed
60+
#define EIDOS_FIND_INTERVAL_USE_BINARY_SEARCH 1
61+
5962
EidosValue *arg_x = p_arguments[0].get();
6063
EidosValue *arg_vec = p_arguments[1].get();
6164
EidosValue *arg_rightmostClosed = p_arguments[2].get();
@@ -72,7 +75,11 @@ EidosValue_SP Eidos_ExecuteFunction_findInterval(const std::vector<EidosValue_SP
7275

7376
bool rightmostClosed = arg_rightmostClosed->LogicalAtIndex(0, nullptr);
7477
bool allInside = arg_allInside->LogicalAtIndex(0, nullptr);
78+
79+
#if !(EIDOS_FIND_INTERVAL_USE_BINARY_SEARCH)
80+
// Used by the old linear-search algorithm
7581
int initial_x_result = allInside ? 0 : -1;
82+
#endif
7683

7784
EidosValue_Int_vector *int_result = (new (gEidosValuePool->AllocateChunk()) EidosValue_Int_vector())->resize_no_initialize(x_count);
7885

@@ -92,6 +99,27 @@ EidosValue_SP Eidos_ExecuteFunction_findInterval(const std::vector<EidosValue_SP
9299
const int64_t *x_data = arg_x->IsSingleton() ? &((EidosValue_Int_singleton *)arg_x)->IntValue_Mutable() : ((EidosValue_Int_vector *)arg_x)->data();
93100

94101
// Find intervals for integer vec, integer x
102+
103+
#if EIDOS_FIND_INTERVAL_USE_BINARY_SEARCH
104+
// binary search - testing indicates this is pretty much always as fast or faster than linear search
105+
for (int x_index = 0; x_index < x_count; ++x_index)
106+
{
107+
int64_t x_value = x_data[x_index];
108+
long x_result;
109+
110+
if (x_value < vec_data[0])
111+
x_result = (allInside ? 0 : -1);
112+
else if (x_value > vec_data[vec_count - 1])
113+
x_result = (allInside ? vec_count - 2 : vec_count - 1);
114+
else if (x_value == vec_data[vec_count - 1])
115+
x_result = ((rightmostClosed || allInside) ? vec_count - 2 : vec_count - 1);
116+
else
117+
x_result = (std::upper_bound(vec_data, vec_data + vec_count, x_value) - vec_data) - 1;
118+
119+
int_result->set_int_no_check(x_result, x_index);
120+
}
121+
#else
122+
// linear search
95123
for (int x_index = 0; x_index < x_count; ++x_index)
96124
{
97125
int64_t x_value = x_data[x_index];
@@ -111,12 +139,34 @@ EidosValue_SP Eidos_ExecuteFunction_findInterval(const std::vector<EidosValue_SP
111139

112140
int_result->set_int_no_check(x_result, x_index);
113141
}
142+
#endif
114143
}
115144
else // (x_type == EidosValueType::kValueFloat)
116145
{
117146
const double *x_data = arg_x->IsSingleton() ? &((EidosValue_Float_singleton *)arg_x)->FloatValue_Mutable() : ((EidosValue_Float_vector *)arg_x)->data();
118147

119148
// Find intervals for integer vec, float x
149+
150+
#if EIDOS_FIND_INTERVAL_USE_BINARY_SEARCH
151+
// binary search - testing indicates this is pretty much always as fast or faster than linear search
152+
for (int x_index = 0; x_index < x_count; ++x_index)
153+
{
154+
double x_value = x_data[x_index];
155+
long x_result;
156+
157+
if (x_value < vec_data[0])
158+
x_result = (allInside ? 0 : -1);
159+
else if (x_value > vec_data[vec_count - 1])
160+
x_result = (allInside ? vec_count - 2 : vec_count - 1);
161+
else if (x_value == vec_data[vec_count - 1])
162+
x_result = ((rightmostClosed || allInside) ? vec_count - 2 : vec_count - 1);
163+
else
164+
x_result = (std::upper_bound(vec_data, vec_data + vec_count, x_value) - vec_data) - 1;
165+
166+
int_result->set_int_no_check(x_result, x_index);
167+
}
168+
#else
169+
// linear search
120170
for (int x_index = 0; x_index < x_count; ++x_index)
121171
{
122172
double x_value = x_data[x_index];
@@ -136,6 +186,7 @@ EidosValue_SP Eidos_ExecuteFunction_findInterval(const std::vector<EidosValue_SP
136186

137187
int_result->set_int_no_check(x_result, x_index);
138188
}
189+
#endif
139190
}
140191
}
141192
else // (vec_type == EidosValueType::kValueFloat))
@@ -154,6 +205,27 @@ EidosValue_SP Eidos_ExecuteFunction_findInterval(const std::vector<EidosValue_SP
154205
const int64_t *x_data = arg_x->IsSingleton() ? &((EidosValue_Int_singleton *)arg_x)->IntValue_Mutable() : ((EidosValue_Int_vector *)arg_x)->data();
155206

156207
// Find intervals for float vec, integer x
208+
209+
#if EIDOS_FIND_INTERVAL_USE_BINARY_SEARCH
210+
// binary search - testing indicates this is pretty much always as fast or faster than linear search
211+
for (int x_index = 0; x_index < x_count; ++x_index)
212+
{
213+
int64_t x_value = x_data[x_index];
214+
long x_result;
215+
216+
if (x_value < vec_data[0])
217+
x_result = (allInside ? 0 : -1);
218+
else if (x_value > vec_data[vec_count - 1])
219+
x_result = (allInside ? vec_count - 2 : vec_count - 1);
220+
else if (x_value == vec_data[vec_count - 1])
221+
x_result = ((rightmostClosed || allInside) ? vec_count - 2 : vec_count - 1);
222+
else
223+
x_result = (std::upper_bound(vec_data, vec_data + vec_count, x_value) - vec_data) - 1;
224+
225+
int_result->set_int_no_check(x_result, x_index);
226+
}
227+
#else
228+
// linear search
157229
for (int x_index = 0; x_index < x_count; ++x_index)
158230
{
159231
int64_t x_value = x_data[x_index];
@@ -173,12 +245,34 @@ EidosValue_SP Eidos_ExecuteFunction_findInterval(const std::vector<EidosValue_SP
173245

174246
int_result->set_int_no_check(x_result, x_index);
175247
}
248+
#endif
176249
}
177250
else // (x_type == EidosValueType::kValueFloat)
178251
{
179252
const double *x_data = arg_x->IsSingleton() ? &((EidosValue_Float_singleton *)arg_x)->FloatValue_Mutable() : ((EidosValue_Float_vector *)arg_x)->data();
180253

181254
// Find intervals for float vec, float x
255+
256+
#if EIDOS_FIND_INTERVAL_USE_BINARY_SEARCH
257+
// binary search - testing indicates this is pretty much always as fast or faster than linear search
258+
for (int x_index = 0; x_index < x_count; ++x_index)
259+
{
260+
double x_value = x_data[x_index];
261+
long x_result;
262+
263+
if (x_value < vec_data[0])
264+
x_result = (allInside ? 0 : -1);
265+
else if (x_value > vec_data[vec_count - 1])
266+
x_result = (allInside ? vec_count - 2 : vec_count - 1);
267+
else if (x_value == vec_data[vec_count - 1])
268+
x_result = ((rightmostClosed || allInside) ? vec_count - 2 : vec_count - 1);
269+
else
270+
x_result = (std::upper_bound(vec_data, vec_data + vec_count, x_value) - vec_data) - 1;
271+
272+
int_result->set_int_no_check(x_result, x_index);
273+
}
274+
#else
275+
// linear search
182276
for (int x_index = 0; x_index < x_count; ++x_index)
183277
{
184278
double x_value = x_data[x_index];
@@ -198,6 +292,7 @@ EidosValue_SP Eidos_ExecuteFunction_findInterval(const std::vector<EidosValue_SP
198292

199293
int_result->set_int_no_check(x_result, x_index);
200294
}
295+
#endif
201296
}
202297
}
203298

0 commit comments

Comments
 (0)