Add kmeans1d benchmark.

#240

bin/kmeans1d.c

@@ -0,0 +1,31 @@
+#include "ccv.h"
+#include <sys/time.h>
+#include <ctype.h>
+#include <3rdparty/dsfmt/dSFMT.h>
+
+static unsigned int get_current_time(void)
+{
+	struct timeval tv;
+	gettimeofday(&tv, NULL);
+	return tv.tv_sec * 1000 + tv.tv_usec / 1000;
+}
+
+int main(int argc, char** argv)
+{
+	dsfmt_t dsfmt;
+	dsfmt_init_gen_rand(&dsfmt, 0);
+	ccv_dense_matrix_t* const a = ccv_dense_matrix_new(1, 1000000, CCV_32F | CCV_C1, 0, 0);
+	int i;
+	for (i = 0; i < 1000000; i++)
+		a->data.f32[i] = dsfmt_genrand_open_close(&dsfmt);
+	int* const clusters = (int*)ccmalloc(sizeof(int) * 1000000);
+	double* const centroids = (double*)ccmalloc(sizeof(double) * 128);
+	unsigned int elapsed_time = get_current_time();
+	ccv_kmeans1d(a, 128, clusters, centroids);
+	elapsed_time = get_current_time() - elapsed_time;
+	printf("elapsed: %ums\n", elapsed_time);
+	ccfree(centroids);
+	ccfree(clusters);
+	ccv_matrix_free(a);
+	return 0;
+}

bin/makefile

@@ -3,7 +3,7 @@ include ../lib/config.mk
 LDFLAGS := -L"../lib" -lccv $(LDFLAGS)
 CFLAGS := -O3 -Wall -I"../lib" $(CFLAGS)
 
-TARGETS = bbffmt msermatch siftmatch bbfcreate bbfdetect scdcreate scddetect swtcreate swtdetect dpmcreate dpmdetect tld icfcreate icfdetect icfoptimize cifar-10 image-net cnnclassify aflw
+TARGETS = bbffmt msermatch siftmatch bbfcreate bbfdetect scdcreate scddetect swtcreate swtdetect dpmcreate dpmdetect tld icfcreate icfdetect icfoptimize cifar-10 image-net cnnclassify aflw kmeans1d
 
 TARGET_SRCS := $(patsubst %,%.c,$(TARGETS))
 

lib/ccv_numeric.c

@@ -1266,21 +1266,65 @@ void ccv_distance_transform(ccv_dense_matrix_t* a, ccv_dense_matrix_t** b, int t
 #undef for_block
 }
 
+__attribute__((__always_inline__))
 inline static double _kmeans1d_cost(double* cumsum, double* cumsum2, int i, int j)
 {
 	if (j < i)
 		return 0;
-	double mu = (cumsum[j + 1] - cumsum[i]) / (j - i + 1);
+	double j_minus_i_plus_1_d = (double)(j - i + 1);
+	double temp = (cumsum[j + 1] - cumsum[i]);
+	double mu = temp / j_minus_i_plus_1_d;
 	double result = cumsum2[j + 1] - cumsum2[i];
-	result += (j - i + 1) * (mu * mu);
-	result -= (2 * mu) * (cumsum[j + 1] - cumsum[i]);
+	result += j_minus_i_plus_1_d * (mu * mu);
+	result -= (2 * mu) * temp;
 	return result;
 }
 
+__attribute__((__always_inline__))
 inline static double _kmeans1d_lookup(double* D, double* cumsum, double* cumsum2, int i, int j)
 {
-	const int col = i < j - 1 ? i : j - 1;
-	return (col >= 0 ? D[col] : 0) + _kmeans1d_cost(cumsum, cumsum2, j, i);
+	const int i_minus_j_plus_1 = i - j + 1;
+	const int col = i_minus_j_plus_1 < 0 ? i : j - 1;
+	double result = (col >= 0 ? D[col] : 0);
+	if (i_minus_j_plus_1 < 1)
+		return result;
+	double i_minus_j_plus_1_d = (double)i_minus_j_plus_1;
+	double temp = (cumsum[i + 1] - cumsum[j]);
+	double mu = temp / i_minus_j_plus_1_d;
+	double result_alt = result + cumsum2[i + 1] - cumsum2[j];
+	result_alt += i_minus_j_plus_1_d * (mu * mu);
+	result_alt -= (2 * mu) * temp;
+	return result_alt;
+}
+
+__attribute__((__always_inline__))
+inline static int _kmeans1d_lookup_compare(double* D, double* cumsum, double* cumsum2, int i, int j1, int j2)
+{
+	// Uses either 2-wide SIMD or instruction-level parallelism.
+	const int i_minus_j1_plus_1 = i - j1 + 1;
+	const int i_minus_j2_plus_1 = i - j2 + 1;
+	const int col1 = i_minus_j1_plus_1 < 0 ? i : j1 - 1;
+	const int col2 = i_minus_j2_plus_1 < 0 ? i : j2 - 1;
+	double result1 = (col1 >= 0 ? D[col1] : 0);
+	double result2 = (col2 >= 0 ? D[col2] : 0);
+
+	double i_minus_j1_plus_1_d = (double)i_minus_j1_plus_1;
+	double i_minus_j2_plus_1_d = (double)i_minus_j2_plus_1;
+	double cumsum_i_1 = cumsum[i + 1];
+	double cumsum2_i_1 = cumsum2[i + 1];
+	double temp1 = (cumsum_i_1 - cumsum[j1]);
+	double temp2 = (cumsum_i_1 - cumsum[j2]);
+	double mu1 = temp1 / i_minus_j1_plus_1_d;
+	double mu2 = temp2 / i_minus_j2_plus_1_d;
+	double result1_alt = result1 + cumsum2_i_1 - cumsum2[j1];
+	double result2_alt = result2 + cumsum2_i_1 - cumsum2[j2];
+	result1_alt += i_minus_j1_plus_1_d * (mu1 * mu1);
+	result2_alt += i_minus_j2_plus_1_d * (mu2 * mu2);
+	result1_alt -= (2 * mu1) * temp1;
+	result2_alt -= (2 * mu2) * temp2;
+	result1 = (i_minus_j1_plus_1 < 1) ? result1 : result1_alt;
+	result2 = (i_minus_j2_plus_1 < 1) ? result2 : result2_alt;
+	return result1 >= result2;
 }
 
 static void _smawk2(int row_start, int row_stride, int row_size, int* cols, int col_size, int* reserved, double* D, double* cumsum, double* cumsum2, int* result)
@@ -1298,7 +1342,7 @@ static void _smawk2(int row_start, int row_stride, int row_size, int* cols, int
 			if (_col_size == 0)
 				break;
 			const int row = row_start + row_stride * (_col_size - 1);
-			if (_kmeans1d_lookup(D, cumsum, cumsum2, row, col) >= _kmeans1d_lookup(D, cumsum, cumsum2, row, _cols[_col_size - 1]))
+			if (_kmeans1d_lookup_compare(D, cumsum, cumsum2, row, col, _cols[_col_size - 1]))
 				break;
 			--_col_size;
 		}
@@ -1352,7 +1396,7 @@ static void _smawk1(int row_start, int row_stride, int row_size, int* cols, int
 			if (_col_size == 0)
 				break;
 			const int row = row_start + row_stride * (_col_size - 1);
-			if (_kmeans1d_lookup(D, cumsum, cumsum2, row, col) >= _kmeans1d_lookup(D, cumsum, cumsum2, row, _cols[_col_size - 1]))
+			if (_kmeans1d_lookup_compare(D, cumsum, cumsum2, row, col, _cols[_col_size - 1]))
 				break;
 			--_col_size;
 		}

lib/nnc/ccv_nnc_easy.h

@@ -287,7 +287,7 @@ static inline void ccv_nnc_tensor_view_get_stride(const ccv_nnc_tensor_view_t* c
 {
 	int x;
 	const int nd = ccv_nnc_tensor_nd(tv->info.dim);
-	const int offset = CCV_NNC_MAX_DIM + 2 - nd;
+	const int offset = ccv_max(CCV_NNC_MAX_DIM + 2 - nd, 0);
 	stride[CCV_NNC_MAX_DIM + 2] = 0;
 	if (CCV_IS_TENSOR_VIEW(tv))
 	{