first entropy calc

TransformationGrid.cu

@@ -34,12 +34,33 @@ float copy_sign(float x,float y){
 __host__
 void TransformationGrid::h_matrix_to_quat(float * const& m,float4 & q){
 	//https://en.wikipedia.org/wiki/Rotation_matrix#Quaternion
-	float t = m[0]*m[0] + m[4]*m[4] + m[8]*m[8];
-	float r = sqrt(1+t);
-	q.w = 0.5*r;
-	q.x = 0.5*(copy_sign(0.5*sqrt(1 + m[0]*m[0] - m[4]*m[4] - m[8]*m[8]),m[7] - m[5]));
-	q.y = 0.5*(copy_sign(0.5*sqrt(1 - m[0]*m[0] + m[4]*m[4] + m[8]*m[8]),m[2] - m[6]));
-	q.z = 0.5*(copy_sign(0.5*sqrt(1 - m[0]*m[0] - m[4]*m[4] + m[8]*m[8]),m[3] - m[1]));
+	float tr = m[0] + m[4] + m[8];
+	float S;
+	if (tr>0) {
+		S = sqrtf(tr + 1)*2;
+		q.w = 0.25*S;
+		q.x = (m[7] - m[5])/S;
+		q.y = (m[2] - m[6])/S;
+		q.z = (m[3] - m[1])/S;
+	} else if ((m[0] > m[4]) && (m[0] > m[8])){
+		S = sqrtf(1 + m[0] - m[4] -m[8])*2;
+		q.w = (m[7] - m[5])/S;
+		q.x = 0.25*S;
+		q.y = (m[1] + m[3])/S;
+		q.z = (m[2] + m[6])/S;
+	} else if (m[4] > m[8]){
+		S = sqrtf(1 + m[4] - m[0] -m[8])*2;
+		q.w = (m[2] - m[6])/S;
+		q.x = (m[1] + m[3])/S;
+		q.y = 0.25*S;
+		q.z = (m[5] + m[7])/S;
+	} else {
+		S = sqrtf(1 + m[8] - m[0] - m[4])*2;
+		q.w = (m[3] - m[1])/S;
+		q.x = (m[2] + m[6])/S;
+		q.y = (m[5] + m[7])/S;
+		q.z = 0.25*S;
+	}
 }
 
 __host__
@@ -839,14 +860,24 @@ struct normalize_chamfer: public thrust::unary_function<float,float>
 struct normalize_envelope: public thrust::unary_function<float,float>
 {
 //	thrust::pair<thrust::device_vector<float>::iterator,thrust::device_vector<float>::iterator> d_minmax;
-    float min;
-    float max;
+    float * d_min;
+    float * d_max;
 	__device__
 	//norm = lambda u: 1/(_max - _min)*u - _min/(_max - _min)
 	//since the chamfer score maps bad solutions to high scores and good solutions to low scores
 	//exchange min am max
 	float operator()(const float& u) const {
-    	return u/(max - min) - min/(max - min);
+    	return u/(d_max[0] - d_min[0]) - d_min[0]/(d_max[0] - d_min[0]);
+	}
+};
+
+
+struct entropy: public thrust::unary_function<float,float>
+{
+	__device__
+	void operator()(float& u) {
+		if (u != 0)
+			u = u*__logf(u);
 	}
 };
 
@@ -936,54 +967,27 @@ void float_min_max(float * d_data,size_t num_data, float * d_min, float * d_max)
 	}
 }
 
-template<size_t block_size>
-__global__
-void shift_and_scale_envelope_score_kernel(float * d_score, float * d_scaled_score, float *d_min, float *d_max, size_t num_scores){
-	__shared__ float reduce[block_size];
-	int num_blocks = num_scores/blockDim.x + 1;
-	if (threadIdx.x == 0 && blockIdx.x == 0){
-		float_min_max<block_size><<<32,block_size>>>(d_score,num_scores,d_min,d_max);
-	}
-	__syncthreads();
-	cudaDeviceSynchronize();
-	int grid_size = blockDim.x*gridDim.x;
-	int g = blockIdx.x*blockDim.x + threadIdx.x;
-	while(g<num_scores){
-		d_scaled_score[g] = d_score[g]/(d_max[0] - d_min[0]) - d_min[0]/(d_max[0] - d_min[0]);
-		g+=grid_size;
-	}
-	__syncthreads();
-	//initiate
-	if (blockIdx.x == 0){
-		reduce[threadIdx.x] = 0;
-	}
-	__syncthreads();
 
-	//t-Design, would usually need weights!
-	int b = blockIdx.x;
-	while (b < num_blocks){
-		int t = b*blockDim.x + threadIdx.x;
-		while (t < num_scores){
-			reduce[threadIdx.x] += d_scaled_score[t];
-			t += blockDim.x;
-		}
-		b += gridDim.x;
-	}
-
-	if (threadIdx.x < 64) reduce[threadIdx.x] += reduce[threadIdx.x + 64];__syncthreads();
-	if (threadIdx.x < 32) reduce[threadIdx.x] += reduce[threadIdx.x + 32];__syncthreads();
-	if (threadIdx.x < 16) reduce[threadIdx.x] += reduce[threadIdx.x + 16];__syncthreads();
-	if (threadIdx.x < 8) reduce[threadIdx.x] += reduce[threadIdx.x + 8];__syncthreads();
-	if (threadIdx.x < 4) reduce[threadIdx.x] += reduce[threadIdx.x + 4];__syncthreads();
-	if (threadIdx.x < 2) reduce[threadIdx.x] += reduce[threadIdx.x + 2];__syncthreads();
-	if (threadIdx.x < 1) reduce[threadIdx.x] += reduce[threadIdx.x + 1];__syncthreads();
-
-	printf("");
-}
 
 __host__
 void TransformationGrid::shift_and_scale_envelope_score(float * d_score, float * d_normalized_score,float *d_min, float *d_max, size_t num_scores){
-	shift_and_scale_envelope_score_kernel<128><<<32,128>>>(d_score,d_normalized_score,d_min,d_max,num_scores);
+//	shift_and_scale_envelope_score_kernel<128><<<32,128>>>(d_score,d_normalized_score,d_min,d_max,num_scores);
+	float_min_max<128><<<32,128>>>(d_score,num_scores,d_min,d_max);
+	normalize_envelope unary_function;
+	unary_function.d_min = d_min;
+	unary_function.d_max = d_max;
+	thrust::transform(thrust::device_pointer_cast(&d_score[0]),thrust::device_pointer_cast(&d_score[num_scores-1]),thrust::device_pointer_cast(&d_normalized_score[0]),unary_function);
+	float * h_sum = (float *) malloc(sizeof(*h_sum));
+//	h_sum[0] = thrust::reduce(thrust::device_pointer_cast(&d_normalized_score[0]),thrust::device_pointer_cast(&d_normalized_score[num_scores-1]));
+	float sum = thrust::reduce(thrust::device_pointer_cast(&d_normalized_score[0]),thrust::device_pointer_cast(&d_normalized_score[num_scores-1]));
+	float * d_sum;
+	cudaMalloc((void **)&d_sum,sizeof(*d_sum));
+	cudaMemcpy(d_sum,h_sum,sizeof(*h_sum),cudaMemcpyHostToDevice);
+	CudaCheckError();
+	thrust::for_each(thrust::device_pointer_cast(&d_normalized_score[0]),thrust::device_pointer_cast(&d_normalized_score[num_scores-1]), thrust::placeholders::_1 /= sum);
+	thrust::for_each(thrust::device_pointer_cast(&d_normalized_score[0]),thrust::device_pointer_cast(&d_normalized_score[num_scores-1]), entropy());
+	float entr = -1.0*thrust::reduce(thrust::device_pointer_cast(&d_normalized_score[0]),thrust::device_pointer_cast(&d_normalized_score[num_scores-1]));
+	printf("Entropy for envelope score: %f\n",entr);
 }
 
 TransformationGrid::~TransformationGrid() {

TransformationGrid.h

@@ -37,6 +37,8 @@
 #include <thrust/execution_policy.h>
 #include <thrust/device_ptr.h>
 #include <thrust/reduce.h>
+#include <thrust/functional.h>
+#include <thrust/for_each.h>
 
 #include <util.h>
 

ccl_test.cu

@@ -580,7 +580,7 @@ BOOST_AUTO_TEST_CASE(envelope_test_simpler)
 		float tau_in = 1.0;
 		std::pair<float3,float3> bounding_box_outer = density.get_bounding_box();
 		std::pair<float3,float3> bounding_box_inner = pixel_coords.get_rotationally_safe_bounding_box();
-		TransformationGrid transformation_grid(bounding_box_outer,bounding_box_inner,binarized_map.h_mid_point,tau_in,"rot_samples/N07_M168_OctaC7.dat");
+		TransformationGrid transformation_grid(bounding_box_outer,bounding_box_inner,binarized_map.h_mid_point,tau_in,"rot_samples/N23_M5880_IcoC7.dat");
 		CudaCheckError();
 		printf("Transformation num: %u\n",transformation_grid.h_num_transformations[0]);
 		printf("Translation num: %u\n",transformation_grid.h_num_translations[0]);