added alignment protocol

AlignmentProtocol.cu

+/*
+ * AlignmentProtocol.cu
+ *
+ *  Created on: Jan 14, 2021
+ *      Author: kkarius
+ */
+
+#include <AlignmentProtocol.h>
+#include <dirent.h>
+#include <stdio.h>
+#include <string.h>
+#include <util.h>
+#include <CRMFtool.h>
+
+#include <IMP/rmf/atom_io.h>
+#include <IMP/rmf.h>
+#include <IMP/atom/Representation.h>
+#include <RMF/FileHandle.h>
+#include <RMF/FileConstHandle.h>
+#include <RMF/infrastructure_macros.h>
+#include <RMF/utility.h>
+
+#include <cuda.h>
+#include <cuda_runtime.h>
+#include <cuda_runtime_api.h>
+
+#pragma pack(1)
+struct M_alpha {
+	float m00;
+	float m01;
+	float m02;
+	float m03;
+	float m11;
+	float m12;
+	float m13;
+	float m22;
+	float m23;
+	float m33;
+	char padding;
+	struct M_alpha& operator+=(const M_alpha& rhs) {
+		m00 += rhs.m00;
+		m01 += rhs.m01;
+		m02 += rhs.m02;
+		m03 += rhs.m03;
+		m11 += rhs.m11;
+		m12 += rhs.m12;
+		m13 += rhs.m13;
+		m22 += rhs.m22;
+		m23 += rhs.m23;
+		m33 += rhs.m33;
+		return *this;}
+};
+
+template <int BLOCKDIM>
+__global__
+void center_and_m_alpha(float * d_particles, float * d_frame0, M_alpha * d_m_alphas, int num_frames, int num_particles){
+	//holds intermediates for center of mass and M_alpha calculations --> size = thread_number*(|M_alpha| + 3*float)
+	extern __shared__ float buffer[];
+	float3 * coms = (float3 *) &buffer[0];
+	M_alpha * malphas = (M_alpha *) &buffer[3*BLOCKDIM];
+
+	//frames
+	int f = 0;
+	//particles
+	int p;
+	//threads
+	int t;
+	//pointer for current frame
+	float * d_frame;
+	//pointer for current particles, plus d for efficiency
+	float4 * p0;
+	float4 * p1;
+	float d;
+	//looping over current frames
+	while (f < num_frames){
+		d_frame = &d_particles[4*f*num_particles];
+		//init coms and malphas to zero
+		t = 0;
+		while (t < BLOCKDIM){
+			coms[t] = {0};
+			malphas[t] = {0};
+			t +=1;
+		}
+		p = threadIdx.x;
+		//looping over particles, calculating center of mass
+		while (p < num_particles){
+			coms[threadIdx.x].x += d_frame[4*p];
+			coms[threadIdx.x].y += d_frame[4*p + 1];
+			coms[threadIdx.x].z += d_frame[4*p + 2];
+			p += BLOCKDIM;
+		}
+		__syncthreads();
+//		if (threadIdx.x < 128){coms[threadIdx.x] += coms[threadIdx.x + 128];}__syncthreads();
+//		if (threadIdx.x < 64){coms[threadIdx.x] += coms[threadIdx.x + 128];}__syncthreads();
+		//here, in effect, BLOCKDIM is forced to be 32. Add statements of the above kind to change this.
+		if (threadIdx.x < 32){
+				coms[threadIdx.x] += coms[threadIdx.x+ 32];
+				coms[threadIdx.x] += coms[threadIdx.x+ 16];
+				coms[threadIdx.x] += coms[threadIdx.x+ 8];
+				coms[threadIdx.x] += coms[threadIdx.x+ 4];
+				coms[threadIdx.x] += coms[threadIdx.x+ 2];
+				coms[threadIdx.x] += coms[threadIdx.x+ 1];
+		}
+		__syncthreads();
+		coms[0] /= __int2float_rn(num_particles);
+		//transform coordinates of frame f to coms
+		p = threadIdx.x;
+		//looping over particles, calculating center of mass
+		while (p < num_particles){
+			d_frame[4*p] -= coms[0].x;
+			d_frame[4*p + 1] -= coms[0].y;
+			d_frame[4*p + 2] -= coms[0].z;
+			p += BLOCKDIM;
+		}
+		__syncthreads();
+		p = threadIdx.x;
+		while (p < num_particles){
+			p0 = (float4 *) &d_frame0[p];
+			p1 = (float4 *) &d_frame[p];
+			d = p1->x * p1->x + p1->y * p1->y + p1->z * p1->z + p0->x * p0->x + p0->y * p0->y
+					+ p0->z * p0->z;
+			malphas[threadIdx.x].m00 = d - 2 * p1->x * p0->x - 2 * p1->y * p0->y - 2 * p1->z * p0->z;
+			malphas[threadIdx.x].m01 = 2 * (p1->y * p0->z - p1->z * p0->y);
+			malphas[threadIdx.x].m02 = 2 * (-p1->x * p0->z + p1->z * p0->x);
+			malphas[threadIdx.x].m03 = 2 * (p1->x * p0->y - p1->y * p0->x);
+			malphas[threadIdx.x].m11 = d - 2 * p1->x * p0->x + 2 * p1->y * p0->y + 2 * p1->z * p0->z;
+			malphas[threadIdx.x].m12 = -2 * (p1->x * p0->y + p1->y * p0->x);
+			malphas[threadIdx.x].m13 = -2 * (p1->x * p0->z + p1->z * p0->x);
+			malphas[threadIdx.x].m22 = d + 2 * p1->x * p0->x - 2 * p1->y * p0->y + 2 * p1->z * p0->z;
+			malphas[threadIdx.x].m23 = -2 * (p1->y * p0->z + p1->z * p0->y);
+			malphas[threadIdx.x].m33 = d + 2 * p1->x * p0->x + 2 * p1->y * p0->y - 2 * p1->z * p0->z;
+			p += BLOCKDIM;
+		}
+		__syncthreads();
+		//reduce malphas
+		//here, in effect, BLOCKDIM is forced to be 32. See remark above.
+		if (threadIdx.x < 32){
+			malphas[threadIdx.x] += malphas[threadIdx.x + 32];
+			malphas[threadIdx.x] += malphas[threadIdx.x + 16];
+			malphas[threadIdx.x] += malphas[threadIdx.x + 8];
+			malphas[threadIdx.x] += malphas[threadIdx.x + 4];
+			malphas[threadIdx.x] += malphas[threadIdx.x + 2];
+			malphas[threadIdx.x] += malphas[threadIdx.x + 1];
+		}
+		d_m_alphas[f] = malphas[0];
+		f += gridDim.x;
+	}
+}
+
+void AlignmentProtocol::_set_reference_frame(void){
+//	if (_rmf_file_paths.size()>0){
+//		RMF::FileConstHandle fh = RMF::open_rmf_file_read_only(_rmf_file_paths[0]);
+//		IMP_NEW(IMP::Model, model, ());
+//		IMP::atom::Hierarchy hierarchy = IMP::rmf::create_hierarchies(fh, model)[0];
+//		IMP::core::GenericHierarchies leaves;
+//			IMP::algebra::Vector3D vec;
+//			double radius;
+//			int f = offset;
+//			for (int file_f = start_frame; file_f < start_frame + amount; file_f++) {
+//				IMP::rmf::load_frame(fh, RMF::FrameID(file_f));
+//				leaves = IMP::core::get_leaves(hierarchy);
+//				IMP::atom::Hierarchy leave_h;
+//				int c = 0;
+//		//		printf("Loading frame %i ...\n", f);
+//				for (auto leave : leaves) {
+//					vec = IMP::core::XYZR(leave).get_coordinates();
+//					radius = IMP::core::XYZR(leave).get_radius();
+//					h_coords_4[f * num_particles + c].x = (float) vec[0];
+//					h_coords_4[f * num_particles + c].y = (float) vec[1];
+//					h_coords_4[f * num_particles + c].z = (float) vec[2];
+//		//			h_coords[4 * (f * ::num_particles + c) + 3] = ::kernel_map[radius];
+//					h_coords_4[f * num_particles + c].w = (float) radius;
+//					c += 1;
+//				}
+//				printf("frame %i loaded into CPU memory ...\n", f);
+//				f++;
+//			}
+//
+//	} else {
+//		printf("Reference frame not set due to the absence of files ...\n");
+//	}
+}
+
+void AlignmentProtocol::_read_dir_entries(void){
+    if (dir_exists(_search_dir)){
+		DIR *di;
+		char *ptr1,*ptr2;
+		int retn;
+		struct dirent *dir;
+		di = opendir(_search_dir.c_str()); //specify the directory name
+		{
+			while ((dir = readdir(di)) != NULL)
+			{
+				ptr1=strtok(dir->d_name,".");
+				ptr2=strtok(NULL,".");
+				if(ptr2!=NULL)
+				{
+					retn = strcmp(ptr2,"rmf");
+					retn = retn && strcmp(ptr2,"rmf3");
+					if(retn==0)
+					{
+//						printf("Found rmf file: %s/%s.%s\n",_search_dir.c_str(),ptr1,ptr2);
+						_rmf_file_paths.push_back(std::string(_search_dir + "/" + std::string(ptr1) + "." + std::string(ptr2)));
+					}
+				}
+			}
+			closedir(di);
+			}
+		}
+    for (std::string _filepath: _rmf_file_paths){
+    	printf("Found rmf file: %s\n",_filepath.c_str());
+    }
+}
+
+void AlignmentProtocol::prepare(){
+	_read_dir_entries();
+}
+
+AlignmentProtocol::~AlignmentProtocol() {
+	// TODO Auto-generated destructor stub
+}
+

AlignmentProtocol.h

+/*
+ * AlignmentProtocol.h
+ *
+ *  Created on: Jan 14, 2021
+ *      Author: kkarius
+ */
+
+#ifndef ALIGNMENTPROTOCOL_H_
+#define ALIGNMENTPROTOCOL_H_
+
+#include<vector>
+#include<string>
+//#include<Particles.h>
+
+class AlignmentProtocol {
+public:
+	AlignmentProtocol(){};
+	AlignmentProtocol(std::string search_dir):_search_dir(search_dir){};
+	void prepare(void);
+	void prepareGPU(int gpu_index);
+	void runOnGPU(int gpu_index);
+	virtual ~AlignmentProtocol();
+private:
+	std::string _search_dir;
+	//the first file will be the reference particle set
+	std::vector<std::string> _rmf_file_paths;
+	void _read_dir_entries(void);
+	void _set_reference_frame(void);
+//	Particles * _reference = new Particles;
+};
+
+#endif /* ALIGNMENTPROTOCOL_H_ */