https://github.com/baharev/ManiSolve
Revision 4df5c5f8c67154d8faf93dafcf5df85ce9899e1a authored by Ali Baharev on 18 November 2019, 10:48:51 UTC, committed by GitHub on 18 November 2019, 10:48:51 UTC
1 parent 8a9951d
Tip revision: 4df5c5f8c67154d8faf93dafcf5df85ce9899e1a authored by Ali Baharev on 18 November 2019, 10:48:51 UTC
Adding a CITATION file
Adding a CITATION file
Tip revision: 4df5c5f
subsampling2.c
// Copyright (C) 2017 University of Vienna
// All rights reserved.
// BSD license.
// Author: Ali Baharev <ali.baharev@gmail.com>
#include <assert.h>
#include <float.h>
#include <math.h>
/* Compile with:
Consider adding -DNDEBUG for further speed:
gcc -c -Ofast -march=native -std=c99 -fPIC -Wall -Wunused-variable -Wextra subsampling2.c
gcc -shared -Ofast subsampling2.o -o subsample2.so
rm subsampling2.o
-----------------------------------------------------------------------
This does not make much sense since there is no allocation in this code:
gcc -c -ggdb3 -fsanitize=address -fno-omit-frame-pointer -O0 -std=c99 -fPIC -Wall -Wextra subsampling.c
gcc -shared -O0 subsampling.o -o subsample.so -Wl,--whole-archive /usr/lib/gcc/x86_64-linux-gnu/5/libasan.a -Wl,--no-whole-archive
rm subsampling.o
No longer required with gcc 4.9.3 and above:
python subsampling.py 2>&1 | misc/va27/asan_symbolize.py
*/
static int max_dmin(const double* d_min,
const int n,
const unsigned char* selected)
{
// Find the max element in the array of the nearest neighbors for the
// not yet selected points.
int i_max = -1;
double d_max = -1.0;
for (int i=0; i<n; ++i) {
if (!selected[i]) {
double d = d_min[i];
if (d > d_max) {
i_max = i;
d_max = d;
}
}
}
assert(i_max != -1);
assert(d_min[i_max] >= 0.0);
return i_max;
}
/*
static double distance(const double* x, int point_dim, int i, int j) {
// Distance in the infinity norm
const double* a = x + i*point_dim;
const double* b = x + j*point_dim;
double dist = 0.0;
for (int k=0; k<point_dim; ++k) {
const double tmp = fabs(a[k] - b[k]);
if (tmp > dist)
dist = tmp;
}
return dist;
}
*/
static double distance(const double* x, int point_dim, int i, int j) {
// Distance in the L1-norm
const double* a = x + i*point_dim;
const double* b = x + j*point_dim;
double dist = 0.0;
for (int k=0; k<point_dim; ++k) {
dist += fabs(a[k] - b[k]);
}
return dist;
}
static void add_points(const double* __restrict x,
unsigned char* __restrict__ selected,
const int n,
const int point_dim,
double* __restrict__ d_min)
{
// Add that non-selected point next whose nearest neighbor among the
// selected points is the most distant.
const int index = max_dmin(d_min, n, selected);
assert(index != -1);
assert(!selected[index]);
selected[index] = 1;
// Update the nearest neighbors. Is the newly selected one the nearest now?
for (int i=0; i<n; ++i) {
if (!selected[i]) { // <-- Can become a bottleneck. Track the selected indices instead then?
const double d = distance(x, point_dim, i, index);
if (d < d_min[i]) {
d_min[i] = d;
}
}
}
}
static void set_d_to_nn(const double* __restrict x, // shape: n x point_dim
const int* __restrict__ idx_selected, // shape: k
const int k,
const int* __restrict__ idx_not_selected, // shape: n-k
const int n, // number of points
const int point_dim,
double* __restrict__ d_min)
{
// Initialize the array of nearest neighbors
for (int ii=0; ii < (n-k); ++ii) {
// for each not selected point
const int i = idx_not_selected[ii];
// find its nearest neighbor among the selected ones
double val_min = DBL_MAX;
for (int jj=0; jj < k; ++jj) {
const int j = idx_selected[jj];
double d = distance(x, point_dim, i, j);
if (d < val_min)
val_min = d;
}
assert(val_min >= 0.0);
assert(val_min < DBL_MAX);
d_min[i] = val_min;
}
}
void downsample2(const double* __restrict__ x, // shape: n x point_dim
const int* __restrict__ idx_selected, // shape: k
int k,
const int* __restrict__ idx_not_selected, // shape: n-k
unsigned char* __restrict__ selected, // shape: n
double* __restrict__ d_min, // shape: n
int n, // number of points
int point_dim,
int n_pts_to_add)
{
// If the calling Python code accidentally passes in a pointer to a view,
// we are reading the wrong data. Try to detect it on a best effort basis.
for (int i=0; i<n*point_dim; ++i)
assert(isfinite(x[i]) &&
"NaN or inf among the point coordinates; x must be contiguous, and not a view.");
// See the calling Python code for the preconditions, they are checked there
set_d_to_nn(x, idx_selected, k, idx_not_selected, n, point_dim, d_min);
for (int i=0; i<n_pts_to_add; ++i) {
add_points(x, selected, n, point_dim, d_min);
}
}
Computing file changes ...