Commit 1cad1604 authored by karius's avatar karius

added rotation grid

parent 0b15ba1b
No preview for this file type
......@@ -254,7 +254,6 @@ void CMrcReader::readData() {
std::fstream::in | std::fstream::binary);
//jump header
x_FileHandle.seekg(1024, std::ios::beg);
x_FileHandle.read((char*) data, bytes);
// double sum = 0;
......
No preview for this file type
File added
No preview for this file type
......@@ -2,6 +2,7 @@
#include <boost/filesystem.hpp>
#include <CudaTest.h>
#include <CRotationGrid.h>
//
//const float feat_R_to_sigma = 0.42466090014400953;
//const float inv_sqrt_2_pi = 0.3989422804014327;
......@@ -33,26 +34,95 @@
// }
//}
int main() {
CRotationGrid rotgrid(
"/data/eclipse-workspace/fitter/cffk-orientation-dc8bb42/data/c48u1.quat",
4);
float * mat_chunk;
#include <stdio.h>
#include <helper_cuda.h>
#define CUDA_ERROR_CHECK
#define CudaSafeCall( err ) __cudaSafeCall( err, __FILE__, __LINE__ )
#define CudaCheckError() __cudaCheckError( __FILE__, __LINE__ )
#define GPU_FRAME_BYTES 4194304
#define THREADS_PER_BLOCK 1024
#define M_PI 3.14159265358979323846
inline void __cudaSafeCall(cudaError err, const char *file, const int line) {
#ifdef CUDA_ERROR_CHECK
if (cudaSuccess != err) {
fprintf( stderr, "cudaSafeCall() failed at %s:%i : %s\n", file, line,
cudaGetErrorString(err));
exit(-1);
}
#endif
return;
}
inline void __cudaCheckError(const char *file, const int line) {
#ifdef CUDA_ERROR_CHECK
cudaError err = cudaGetLastError();
if (cudaSuccess != err) {
fprintf( stderr, "cudaCheckError() failed at %s:%i : %s\n", file, line,
cudaGetErrorString(err));
exit(-1);
}
// More careful checking. However, this will affect performance.
// Comment away if needed.
err = cudaDeviceSynchronize();
if (cudaSuccess != err) {
fprintf( stderr, "cudaCheckError() with sync failed at %s:%i : %s\n",
file, line, cudaGetErrorString(err));
exit(-1);
}
#endif
return;
}
texture<float, 1, cudaReadModeElementType> texInput;
__global__
void copyKernel(float*output, int n) {
for (int i = 0; i < n; i++) {
output[i] = tex1Dfetch(texInput, i);
}
}
int main(int argc, char*argv[]) {
const int WIDTH = 8;
float* hInput = (float*) malloc(sizeof(float) * WIDTH);
float*hOutput = (float*) malloc(sizeof(float) * WIDTH);
for (int i = 0; i < WIDTH; i++) {
hInput[i] = (float) i;
}
float* dInput = NULL, *dOutput = NULL;
size_t offset = 0;
while (!rotgrid) {
mat_chunk = *rotgrid;
for (int i = 0; i < 4 * 9; i++) {
printf("%f\t", mat_chunk[i]);
if ((i + 1) % 3 == 0) {
printf("\n");
}
if ((i + 1) % 9 == 0) {
printf("=============================\n");
}
}
printf("CHUNK===============================\n");
rotgrid++;
texInput.addressMode[0] = cudaAddressModeBorder;
texInput.addressMode[1] = cudaAddressModeBorder;
texInput.filterMode = cudaFilterModePoint;
texInput.normalized = false;
checkCudaErrors(cudaMalloc((void** )&dInput, sizeof(float) * WIDTH));
checkCudaErrors(cudaMalloc((void** )&dOutput, sizeof(float) * WIDTH));
cudaMemcpy(dInput, hInput, sizeof(float) * WIDTH, cudaMemcpyHostToDevice);
cudaBindTexture(&offset, texInput, dInput, sizeof(float) * WIDTH);
copyKernel<<<1,1>>>(dOutput, WIDTH);
cudaMemcpy(hOutput, dOutput, sizeof(float) * WIDTH, cudaMemcpyDeviceToHost);
printf("\nInput = ");
for (int i = 0; i < WIDTH; i++) {
printf("%f\t", hInput[i]);
}
printf("\nOutput = ");
for (int i = 0; i < WIDTH; i++) {
printf("%f\t", hOutput[i]);
}
return 0;
......
......@@ -13,8 +13,11 @@ CUDA_LINKS = -L/usr/local/cuda-10.0/lib64/ -lcusolver -lcudart -lcublas
LFLAGS=-L. $(IMP_LINKS) $(BOOST_LINKS) $(CUDA_LINKS)
all:
$(NVCC) $(CFLAGS) $(LFLAGS) testareal.cu -o testareal
all: CRotationGrid.o
$(NVCC) $(CFLAGS) $(LFLAGS) testareal.cu -o testareal CRotationGrid.o
CRotationGrid.o: CRotationGrid.cpp
$(NVCC) -c -o $@ $< $(CFLAGS)
This diff is collapsed.
This diff is collapsed.
This diff is collapsed.
This diff is collapsed.
This diff is collapsed.
This diff is collapsed.
This diff is collapsed.
This diff is collapsed.
This diff is collapsed.
This diff is collapsed.
This diff is collapsed.
This diff is collapsed.
This diff is collapsed.
This diff is collapsed.
This diff is collapsed.
This diff is collapsed.
This diff is collapsed.
This diff is collapsed.
This diff is collapsed.
This diff is collapsed.
This diff is collapsed.
This diff is collapsed.
This diff is collapsed.
This diff is collapsed.
This diff is collapsed.
This diff is collapsed.
#/usr/bin/python
import RMF,IMP,IMP.rmf
fh = RMF.open_rmf_file_read_only("./concat.rmf3")
m = IMP.Model()
h = IMP.rmf.create_hierarchies(fh,m)[0]
ofh=RMF.create_rmf_file('./py_concat_n.rmf')
IMP.rmf.add_hierarchy(ofh,h)
rbs = [e for e in IMP.atom.get_leaves(h) if IMP.core.RigidBody_get_is_setup(e)]
for rb in rbs:
IMP.core.transform(IMP.core.RigidBody(rb),IMP.algebra.Transformation3D(IMP.algebra.Vector3D(-30,-30,-30)))
m.update()
IMP.rmf.save_frame(ofh)
for rb in rbs:
IMP.core.transform(IMP.core.RigidBody(rb),IMP.algebra.Transformation3D(IMP.algebra.Vector3D(-30,-30,-30)))
m.update()
IMP.rmf.save_frame(ofh)
for rb in rbs:
IMP.core.transform(IMP.core.RigidBody(rb),IMP.algebra.Transformation3D(IMP.algebra.Vector3D(-30,-30,-30)))
m.update()
IMP.rmf.save_frame(ofh)
del ofh
No preview for this file type
......@@ -18,6 +18,8 @@
#include <assert.h>
#include <cmath>
#include <cuda_util_include/cutil_math.h>
#include <CRotationGrid.h>
#define CUDA_ERROR_CHECK
#define CudaSafeCall( err ) __cudaSafeCall( err, __FILE__, __LINE__ )
......@@ -38,6 +40,27 @@ inline void __cudaSafeCall(cudaError err, const char *file, const int line) {
return;
}
inline void __cudaCheckError(const char *file, const int line) {
#ifdef CUDA_ERROR_CHECK
cudaError err = cudaGetLastError();
if (cudaSuccess != err) {
fprintf( stderr, "cudaCheckError() failed at %s:%i : %s\n", file, line,
cudaGetErrorString(err));
exit(-1);
}
// More careful checking. However, this will affect performance.
// Comment away if needed.
err = cudaDeviceSynchronize();
if (cudaSuccess != err) {
fprintf( stderr, "cudaCheckError() with sync failed at %s:%i : %s\n",
file, line, cudaGetErrorString(err));
exit(-1);
}
#endif
return;
}
size_t num_frames;
size_t num_particles;
......@@ -167,79 +190,6 @@ float * getRotationMatrices(std::string quat_file,size_t & num_rotations) {
return matrices;
}
float * iterRotationMatrices(std::ifstream fs,float * matrices,size_t iter_num) {
;
;
fs.open(quat_file, std::fstream::in);
std::string line = "";
bool found_start = false;
std::vector<std::string> strs;
size_t n_rotations;
while (!fs.eof()) {
std::getline(fs, line);
// std::cout << line << std::endl;
if (line.find("format") != std::string::npos) {
//skip parameter lines
std::getline(fs, line);
std::getline(fs, line);
boost::split(strs, line, boost::is_any_of(" \t"),
boost::token_compress_on);
// std::cout << strs[0] << std::endl;
n_rotations = std::stoi(strs[0]);
num_rotations = n_rotations;
found_start = true;
break;
}
}
matrices = (float *) malloc(n_rotations * 9 * sizeof(float));
float * matrix = (float *) malloc(9 * sizeof(float));
float * quat = (float*) malloc(4 * sizeof(float));
if (found_start) {
int rot_offset = 0;
while (!fs.eof()) {
std::getline(fs, line);
boost::trim(line);
if (line.size() != 0) {
boost::split(strs, line, boost::is_any_of(" \t"),
boost::token_compress_on);
if (strs.size() >= 4) {
for (int i = 0; i < 4; i++) {
quat[i] = std::stof(strs[i]);
}
}
matrix = getRotationMatrix(quat);
}
for (int i = 0; i < 9; i++) {
matrices[rot_offset + i] = matrix[i];
}
rot_offset += 9;
}
}
return matrices;
}
inline void __cudaCheckError(const char *file, const int line) {
#ifdef CUDA_ERROR_CHECK
cudaError err = cudaGetLastError();
if (cudaSuccess != err) {
fprintf( stderr, "cudaCheckError() failed at %s:%i : %s\n", file, line,
cudaGetErrorString(err));
exit(-1);
}
// More careful checking. However, this will affect performance.
// Comment away if needed.
err = cudaDeviceSynchronize();
if (cudaSuccess != err) {
fprintf( stderr, "cudaCheckError() with sync failed at %s:%i : %s\n",
file, line, cudaGetErrorString(err));
exit(-1);
}
#endif
return;
}
//IMP::atom::Hierarchy get_hierarchy(std::string rmf_path, IMP::Model * m) {
// //figure out how to get rid of file handle
......@@ -742,10 +692,10 @@ bool test_reduced_coords(float4 * d_reduced_coords, size_t num_frames) {
return ret;
}
void calc_rmsds(float4 * d_coords_4, float * d_rmsds,size_t num_particles, size_t num_frames,size_t threads_per_block,size_t cycle_vol){
void calc_rmsds(float4 * d_coords_4, float * d_rmsds,size_t num_particles, size_t num_frames,size_t threads_per_block){
size_t rmsd_iterations = num_particles / threads_per_block + 1;
for (int i = 0; i < rmsd_iterations; i++) {
calc_rmsd<<<cycle_vol, threads_per_block,sizeof(float)*threads_per_block>>>(d_coords_4,num_particles,i*threads_per_block,d_rmsds);
calc_rmsd<<<num_frames, threads_per_block,sizeof(float)*threads_per_block>>>(d_coords_4,num_particles,i*threads_per_block,d_rmsds);
}
calc_rmsd_sqrt<<<1,num_frames>>>(d_rmsds,num_frames);
}
......@@ -1113,68 +1063,58 @@ int main() {
cudaDeviceSynchronize();
printf("Copying done ...\n");
//<----------------- RMSD TEST ----------------------->
size_t n_rots;
float * h_test_rotations = getRotationMatrices("/data/eclipse-workspace/fitter/cffk-orientation-dc8bb42/data/c48n9.quat",n_rots);
float * rmsds = (float *) malloc(n_rots*sizeof(float));
float * d_test_rotations;
cudaMalloc(&d_test_rotations,sizeof(float)*9*n_rots);
cudaMemcpy(d_test_rotations,h_test_rotations,sizeof(float)*9*n_rots,cudaMemcpyHostToDevice);
//copy the one frame needed to second entry
cudaMemcpy(d_coords+4*num_particles,d_transformed_coords+4*num_particles,4*num_particles*sizeof(float),cudaMemcpyDeviceToDevice);
for (int r = 0;r<cycle_vol-1;r++){
for (int c = 0; c < num_particles; c++){
cublas_status = cublasSgemv(handle, CUBLAS_OP_N, 3, 3, &alpha,
d_test_rotations + 9 * r, 3,
d_coords + 4*(num_particles + c), 1, &beta,
d_coords + 4 *((r +1) * num_particles + c),
1);
}
}
print_d_array(d_coords,5,5,4*num_particles,4*num_frames*num_particles);
if (false){
//copy the one frame needed to second entry
cudaMemcpy(d_coords+4*num_particles,d_transformed_coords+4*num_particles,4*num_particles*sizeof(float),cudaMemcpyDeviceToDevice);
size_t rotgrid_chunksize = 36;
CRotationGrid rotgrid("/data/eclipse-workspace/fitter/cffk-orientation-dc8bb42/data/c48u4749.quat",rotgrid_chunksize);
float * d_rotsample;
cudaMalloc(&d_rotsample,rotgrid_chunksize*9*sizeof(float));
float * d_rmsds;
cudaMalloc(&d_rmsds, sizeof(float) * num_frames);
float * h_rmsds = (float *) malloc(sizeof(float)*rotgrid_chunksize);
bool test_success = true;
size_t ch = 0;
while(!rotgrid){
// rotsample = *rotgrid;
cudaMemcpy(d_rotsample,*rotgrid,rotgrid_chunksize*9*sizeof(float),cudaMemcpyHostToDevice);
for (int r = 0;r<rotgrid_chunksize;r++){
for (int c = 0; c < num_particles; c++){
//first rotation of first chunk is identity
if (r==0 && ch==0)
continue;
cublas_status = cublasSgemv(handle, CUBLAS_OP_N, 3, 3, &alpha,
d_rotsample + 9 * r, 3,
d_coords + 4*(num_particles + c), 1, &beta,
d_coords + 4 *((r +2) * num_particles + c),
1);
}
}
//copy rmsds to host to test them
calc_rmsds(d_coords_4,d_rmsds,num_particles,rotgrid_chunksize,THREADS_PER_BLOCK);
cudaMemcpy(h_rmsds,d_rmsds,sizeof(float)*rotgrid_chunksize,cudaMemcpyDeviceToHost);
for (int j=1;j<rotgrid_chunksize;j++){
test_success &= h_rmsds[0] < h_rmsds[j];
}
float * d_rmsds;
cudaMalloc(&d_rmsds, sizeof(float) * num_frames);
calc_rmsds(d_coords_4,d_rmsds,num_particles,cycle_vol,THREADS_PER_BLOCK,cycle_vol);
printf("Testing orientations %i to %i: %s ...\n",ch*rotgrid_chunksize,(ch+1)*rotgrid_chunksize,test_success ? "success":"failure" );
ch++;
rotgrid++;
}
// print_d_array(d_coords,5,5,4*num_particles,4*num_frames*num_particles);
}
float * h_rmsds = (float *) malloc(sizeof(float)*num_frames);
cudaMemcpy(h_rmsds,d_rmsds,sizeof(float)*cycle_vol,cudaMemcpyDeviceToHost);
// cudaDeviceSynchronize();
CudaCheckError();
//
printf("rmsd[0]:%f\n",h_rmsds[0]);
printf("rmsd[1]:%f\n",h_rmsds[1]);
printf("rmsd[2]:%f\n",h_rmsds[2]);
printf("rmsd[3]:%f\n",h_rmsds[3]);
printf("rmsd[4]:%f\n",h_rmsds[4]);
cublas_status = cublasDestroy(handle);
printf("Status cublas context destruction: %s \n",
cublas_status ? "[failure] ..." : "[success] ..." );
cudaDeviceSynchronize();
CudaCheckError();
// RMF::FileHandle ofh = RMF::create_rmf_file(
// "/data/eclipse-workspace/fitter/test/concat_out.rmf3");
// IMP::rmf::add_hierarchy(ofh, hierarchy);
// IMP::core::GenericHierarchies leaves = IMP::core::get_leaves(hierarchy);
// IMP::algebra::Vector3D vec;
// for (int f = 0; f < cycle_vol; f++) {
// size_t c = 0;
// for (auto leave : leaves) {
// vec[0] = frames[3 * (f * num_particles + c) + 0];
// vec[1] = frames[3 * (f * num_particles + c) + 1];
// vec[2] = frames[3 * (f * num_particles + c) + 2];
// IMP::core::XYZR(leave).set_coordinates(vec);
// c++;
// }
//// printf("Writing frame number: %i ...\n", f);
// IMP::rmf::save_frame(ofh);
// ofh.flush();
// }
return 0;
}
Markdown is supported
0%
or
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment