#include <cmath>
#include <string>
#include <sstream>
+#include <cfloat>
#include <gromacs/trajectoryanalysis.h>
#include <gromacs/utility/smalloc.h>
using gmx::RVec;
+double F (double aix, double aiy, double aiz, double bix_plus_x, double biy_plus_y, double biz_plus_z, double A, double B, double C) {
+ return sqrt( pow(aix + biy_plus_y * (cos(A) * sin(C) - cos(C) * sin(A) * sin(B)) - biz_plus_z * (sin(A) * sin(C) + cos(A) * cos(C) * sin(B)) - cos(B) * cos(C) * bix_plus_x, 2) +
+ pow(aiy - biy_plus_y * (cos(A) * cos(C) + sin(A) * sin(B) * sin(C)) + biz_plus_z * (cos(C) * sin(A) - cos(A) * sin(B) * sin(C)) - cos(B) * sin(C) * bix_plus_x, 2) +
+ pow(aiz + sin(B) * bix_plus_x - cos(A) * cos(B) * biz_plus_z - cos(B) * sin(A) * biy_plus_y, 2) );
+}
+
+void searchF0xyzabc (double &F, double &Fx, double &Fy, double &Fz, double &Fa, double &Fb, double &Fc, double aix, double aiy, double aiz, double bix_plus_x, double biy_plus_y, double biz_plus_z, double A, double B, double C) {
+ double t01, t02, t03, t04, t05, t06, t07;
+ t01 = pow(aix + biy_plus_y * (cos(A) * sin(C) - cos(C) * sin(A) * sin(B)) - biz_plus_z * (sin(A) * sin(C) + cos(A) * cos(C) * sin(B)) - cos(B) * cos(C) * bix_plus_x, 2);
+ t02 = pow(aiy - biy_plus_y * (cos(A) * cos(C) + sin(A) * sin(B) * sin(C)) + biz_plus_z * (cos(C) * sin(A) - cos(A) * sin(B) * sin(C)) - cos(B) * sin(C) * bix_plus_x, 2);
+ t03 = pow(aiz + sin(B) * bix_plus_x - cos(A) * cos(B) * biz_plus_z - cos(B) * sin(A) * biy_plus_y, 2);
+ t04 = sqrt(t01 + t02 + t03);
+ t05 = (aix + biy_plus_y * (cos(A) * sin(C) - cos(C) * sin(A) * sin(B)) - biz_plus_z * (sin(A) * sin(C) + cos(A) * cos(C) * sin(B)) - cos(B) * cos(C) * bix_plus_x);
+ t06 = (aiz + sin(B) * bix_plus_x - cos(A) * cos(B) * biz_plus_z - cos(B) * sin(A) * biy_plus_y);
+ t07 = (aiy - biy_plus_y * (cos(A) * cos(C) + sin(A) * sin(B) * sin(C)) + biz_plus_z * (cos(C) * sin(A) - cos(A) * sin(B) * sin(C)) - cos(B) * sin(C) * bix_plus_x);
+ F += t04;
+ Fx += -(2 * cos(B) * cos(C) * t05 - 2 * sin(B) * t06 + 2 * cos(B) * sin(C) * t07) / (2 * t04);
+ Fy += -(2 * (cos(A) * cos(C) + sin(A) * sin(B) * sin(C)) * t07 - 2 * (cos(A) * sin(C) - cos(C) * sin(A) * sin(B)) * t05 + 2 * cos(B) * sin(A) * t06) / (2 * t04);
+ Fz += -(2 * (sin(A) * sin(C) + cos(A) * cos(C) * sin(B)) * t05 - 2 * (cos(C) * sin(A) - cos(A) * sin(B) * sin(C)) * t07 + 2 * cos(A) * cos(B) * t06) / (2 * t04);
+ Fa += -(2 * (cos(A) * cos(B) * biy_plus_y - cos(B) * sin(A) * biz_plus_z) * t06 -
+ 2 * (biy_plus_y * (cos(C) * sin(A) - cos(A) * sin(B) * sin(C)) + biz_plus_z * (cos(A) * cos(C) + sin(A) * sin(B) * sin(C))) * t07 +
+ 2 * (biy_plus_y * (sin(A) * sin(C) + cos(A) * cos(C) * sin(B)) + biz_plus_z * (cos(A) * sin(C) - cos(C) * sin(A) * sin(B))) * t05) / (2 * t04);
+ Fb += -(2 * (cos(A) * cos(B) * sin(C) * biz_plus_z - sin(B) * sin(C) * bix_plus_x + cos(B) * sin(A) * sin(C) * biy_plus_y) * t07 +
+ 2 * (cos(A) * cos(B) * cos(C) * biz_plus_z - cos(C) * sin(B) * bix_plus_x + cos(B) * cos(C) * sin(A) * biy_plus_y) * t05 -
+ 2 * (cos(B) * bix_plus_x + sin(A) * sin(B) * biy_plus_y + cos(A) * sin(B) * biz_plus_z) * t06) / (2 * t04);
+ Fc += (2 * (biy_plus_y * (cos(A) * cos(C) + sin(A) * sin(B) * sin(C)) - biz_plus_z * (cos(C) * sin(A) - cos(A) * sin(B) * sin(C)) + cos(B) * sin(C) * bix_plus_x) * t05 -
+ 2 * (biz_plus_z * (sin(A) * sin(C) + cos(A) * cos(C) * sin(B)) - biy_plus_y * (cos(A) * sin(C) - cos(C) * sin(A) * sin(B)) + cos(B) * cos(C) * bix_plus_x) * t07) / (2 * t04);
+}
+
+void ApplyFit (std::vector< RVec > &b, double x, double y, double z, double A, double B, double C) {
+ double t0 = 0, t1 = 0, t2 = 0;
+ for (unsigned int i = 0; i < b.size(); i++) {
+ t0 = static_cast< double >(b[i][0]);
+ t1 = static_cast< double >(b[i][1]);
+ t2 = static_cast< double >(b[i][2]);
+ b[i][0] = static_cast< float >((t2 + z) * (sin(A) * sin(C) + cos(A) * cos(C) * sin(B)) - (t1 + y) * (cos(A) * sin(C) - cos(C) * sin(A) * sin(B)) + cos(B) * cos(C) * (t0 + x));
+ b[i][1] = static_cast< float >((t1 + y) * (cos(A) * cos(C) + sin(A) * sin(B) * sin(C)) - (t2 + z) * (cos(C) * sin(A) - cos(A) * sin(B) * sin(C)) + cos(B) * sin(C) * (t0 + x));
+ b[i][2] = static_cast< float >(cos(A) * cos(B) * (t2 + z) - sin(B) * (t0 + x) + cos(B) * sin(A) * (t1 + y));
+ }
+}
+
+void CalcMid (std::vector< RVec > a, std::vector< RVec > b, RVec &midA, RVec &midB, std::vector< std::pair< unsigned int, unsigned int > > pairs) {
+ midA[0] = 0;
+ midA[1] = 0;
+ midA[2] = 0;
+
+ midB[0] = 0;
+ midB[1] = 0;
+ midB[2] = 0;
+
+ for (unsigned int i = 0; i < pairs.size(); i++) {
+ //midA += a[pairs[i].first];
+ //midB += b[pairs[i].second];
+ rvec_inc(midA, a[pairs[i].first]);
+ rvec_inc(midB, b[pairs[i].second]);
+ }
+ midA[0] /= pairs.size();
+ midA[1] /= pairs.size();
+ midA[2] /= pairs.size();
+
+ midB[0] /= pairs.size();
+ midB[1] /= pairs.size();
+ midB[2] /= pairs.size();
+}
+
+void MyFitNew (std::vector< RVec > a, std::vector< RVec > &b, std::vector< std::pair< unsigned int, unsigned int > > pairs, double FtCnst) {
+ double f1 = 0, f2 = 0, fx = 0, fy = 0, fz = 0, fa = 0, fb = 0, fc = 0, l = 1;
+ RVec ma, mb;
+ CalcMid(a, b, ma, mb, pairs);
+ //ma -= mb;
+ rvec_dec(ma, mb);
+ double x = static_cast< double >(ma[0]), y = static_cast< double >(ma[1]), z = static_cast< double >(ma[2]), A = 0, B = 0, C = 0;
+ for (unsigned int i = 0; i < pairs.size(); i++) {
+ f1 += F(static_cast< double >(a[pairs[i].first][0]), static_cast< double >(a[pairs[i].first][1]), static_cast< double >(a[pairs[i].first][2]),
+ static_cast< double >(b[pairs[i].second][0]) + x, static_cast< double >(b[pairs[i].second][1]) + y, static_cast< double >(b[pairs[i].second][2]) + z, A, B, C);
+ }
+ if (FtCnst == 0) {
+ FtCnst = 0.00001;
+ }
+ if (f1 == 0) {
+ return;
+ } else {
+ while (f1 - f2 < 0 || f1 - f2 > FtCnst) {
+ f1 = 0; fx = 0; fy = 0; fz = 0; fa = 0; fb = 0; fc = 0; l = 1;
+ for (unsigned int i = 0; i < pairs.size(); i++) {
+ searchF0xyzabc(f1, fx, fy, fz, fa, fb, fc,
+ static_cast< double >(a[pairs[i].first][0]), static_cast< double >(a[pairs[i].first][1]), static_cast< double >(a[pairs[i].first][2]),
+ static_cast< double >(b[pairs[i].second][0]) + x, static_cast< double >(b[pairs[i].second][1]) + y,
+ static_cast< double >(b[pairs[i].second][2]) + z, A, B, C);
+ }
+ while (f2 != f1) {
+ f2 = 0;
+ for (unsigned int i = 0; i < pairs.size(); i++) {
+ f2 += F(static_cast< double >(a[pairs[i].first][0]), static_cast< double >(a[pairs[i].first][1]), static_cast< double >(a[pairs[i].first][2]),
+ static_cast< double >(b[pairs[i].second][0]) + (x - l * fx), static_cast< double >(b[pairs[i].second][1]) + (y - l * fy),
+ static_cast< double >(b[pairs[i].second][2]) + (z - l * fz), A - l * fa, B - l * fb, C - l * fc);
+ }
+ if (f2 < f1) {
+ x -= l * fx; y -= l * fy; z -= l * fz; A -= l * fa; B -= l * fb; C -= l * fc;
+ fx = 0; fy = 0; fz = 0; fa = 0; fb = 0; fc = 0;
+ break;
+ } else {
+ l *= 0.50;
+ if (l == DBL_MIN) { //DBL_TRUE_MIN
+ break;
+ }
+ }
+ }
+ }
+ ApplyFit(b, x, y, z, A, B, C);
+ }
+}
+
+//const int tau_jump = 10;
+
const std::vector< std::vector< std::vector< double > > > read_correlation_matrix_file(const char* file_name, unsigned long size, unsigned int length)
{
FILE *file;
std::vector< std::vector< double > > b;
b.resize(size, c);
crrlts.resize(length + 1, b);
- char a[100];
+ //char a[100];
std::cout << "reading correlations from a file:\n";
for (unsigned int i = 0; i < length + 1; i++) {
int t0 = 0, t1 = std::fscanf(file, "%d\n", &t0);
}
template< typename T >
-unsigned int posSrch(std::vector< T > a, T b) {
+unsigned long posSrch(std::vector< T > a, T b) {
for (unsigned i01 = 0; i01 < a.size(); i01++) {
if (a[i01] == b) {
return i01;
void make_rout_file(double crl_border, std::vector< int > indx, std::vector< std::string > resNames,
std::vector< std::vector< std::vector< std::pair< unsigned int, unsigned int > > > > rout,
- const char* file_name01, const unsigned int tauD, const unsigned int window_size)
+ const char* file_name01, const unsigned int tauD, const unsigned int window_size, int t_jump)
{
FILE *file01;
file01 = std::fopen(file_name01, "w+");
}
}
if (i01 == f) {
- std::fprintf(file01, "\n Starting time point = %d | correlations >= %0.2f | tau = %d | window = %d\n\n", i01, crl_border, tauD, window_size);
+ std::fprintf(file01, "\n Starting time point = %d | correlations >= %0.2f | tau = %d | window = %d\n\n", i01 * window_size / t_jump, crl_border, tauD, window_size);
} else {
- std::fprintf(file01, "\n Starting time point = [%d ; %d] | correlations >= %0.2f | tau = %d | window = %d\n\n", i01, f, crl_border, tauD, window_size);
+ std::fprintf(file01, "\n Starting time point = [%d ; %d] | correlations >= %0.2f | tau = %d | window = %d\n\n", i01 * window_size / t_jump, f * window_size / t_jump, crl_border, tauD, window_size);
}
for (unsigned int i02 = 0; i02 < rout[i01].size(); i02++) {
for (unsigned int i03 = 0; i03 < rout[i01][i02].size(); i03++) {
table.resize(0);
for (unsigned int i02 = 0; i02 < rout[i01].size(); i02++) {
for (unsigned int i03 = 0; i03 < rout[i01][i02].size(); i03++) {
- bool flag = true;
+ bool flag1 = true, flag2 = true;
for (unsigned int i04 = 0; i04 < table.size(); i04++) {
if (std::get<0>(table[i04]) == indx[rout[i01][i02][i03].first]) {
std::get<1>(table[i04])++;
std::get<2>(table[i04]).push_back(indx[rout[i01][i02][i03].second]);
- flag = false;
- } else if (std::get<0>(table[i04]) == indx[rout[i01][i02][i03].second]) {
+ flag1 = false;
+ }
+ if (std::get<0>(table[i04]) == indx[rout[i01][i02][i03].second]) {
std::get<1>(table[i04])++;
std::get<2>(table[i04]).push_back(indx[rout[i01][i02][i03].first]);
- flag = false;
+ flag2 = false;
}
}
- if (flag) {
+ if (flag1) {
std::vector< int > a;
a.resize(0);
a.push_back(indx[rout[i01][i02][i03].second]);
table.push_back(std::make_tuple(indx[rout[i01][i02][i03].first], 1, a));
+ }
+ if (flag2) {
+ std::vector< int > a;
a.resize(0);
a.push_back(indx[rout[i01][i02][i03].first]);
table.push_back(std::make_tuple(indx[rout[i01][i02][i03].second], 1, a));
}
}
-void correlation_evaluation(const std::vector< RVec > ref, const std::vector< std::vector< RVec > > trajectory, const unsigned int tauE, const unsigned int window_size, const char* file_name) {
+float minDistToSubset(const std::vector< RVec > set, std::vector < unsigned int > subSet, int pos) {
+ float a = 9000;
+ for (unsigned int j = 0; j < subSet.size(); j++) {
+ if (a > gmx::norm((set[pos]-set[j]).as_vec())) {
+ a = gmx::norm((set[pos]-set[j]).as_vec());
+ }
+ }
+ return a;
+}
+
+void correlation_evaluation(const std::vector< RVec > ref,
+ const std::vector< std::vector< RVec > > trajectory,
+ const unsigned int tauE, const unsigned int window_size, const char* file_name, int t_jump,
+ std::vector< std::vector < std::vector < unsigned int > > > sls) {
std::vector< std::vector< std::vector< double > > > crl;
std::vector< std::vector< RVec > > trj;
- for (unsigned int istart = 0; istart < trajectory.size() - window_size - tauE; istart++) {
+ for (unsigned int istart = 0; istart < trajectory.size() - window_size - tauE; istart += window_size / t_jump) {
trj.clear();
trj.resize(window_size + tauE);
for (unsigned int i = 0; i < window_size + tauE; i++) {
trj[i].push_back(trajectory[i + istart][j]);
}
}
+
+ /*
+ *
+ * fitting block
+ * we add spare atoms to existing domains based on proximity
+ *
+ */
+
+ std::vector< std::vector< std::pair< unsigned int, unsigned int > > > prs;
+ int prsCount = 0;
+ float prsTemp = 9000;
+ for (unsigned int i0 = 0; i0 < ref.size(); i0++) {
+ for (unsigned int i = 0; i < sls[istart / (window_size / t_jump)].size(); i++) {
+ if (prsTemp > minDistToSubset(ref, sls[istart / (window_size / t_jump)][i], i0)) {
+ prsTemp = minDistToSubset(ref, sls[istart / (window_size / t_jump)][i], i0);
+ prsCount = i;
+ }
+ if (prsTemp < 0.1) { // to avoid 0. error
+ prsCount = -1;
+ break;
+ }
+ }
+ if (prsCount != -1) {
+ sls[istart / (window_size / t_jump)][prsCount].push_back(i0);
+ }
+ }
+
+ prs.resize(sls[istart / (window_size / t_jump)].size());
+ for (unsigned int i = 0; i < sls[istart / (window_size / t_jump)].size(); i++) {
+ for (unsigned int j = 0; j < sls[istart / (window_size / t_jump)][i].size(); j++) {
+ prs[i].push_back(std::make_pair(sls[istart / (window_size / t_jump)][i][j], sls[istart / (window_size / t_jump)][i][j]));
+ }
+ }
+
+ std::vector< std::vector< std::vector< RVec > > > trjTemp;
+ trjTemp.resize(prs.size());
+ for (unsigned int i = 0; i < prs.size(); i++) {
+ trjTemp[i] = trj;
+ for (unsigned int j = 0; j < trjTemp[i].size(); j++) {
+ MyFitNew(ref, trjTemp[i][j], prs[i], 0);
+ for (unsigned int k = 0; k < prs[i].size(); k++) {
+ trj[j][prs[i][k].first] = trjTemp[i][j][prs[i][k].first];
+ }
+ }
+ }
+
+ /*
+ *
+ * fitting block
+ *
+ */
+
+
crl.resize(tauE + 1);
for (unsigned int i = 0; i < crl.size(); i++) {
crl[i].resize(trajectory.front().size());
std::vector< std::vector< double > > a, b, c;
std::vector< double > d;
d.resize(0);
+ a.resize(0);
+ b.resize(0);
+ c.resize(0);
for (unsigned int j = 0; j < trajectory.front().size(); j++) {
a.push_back(d);
b.push_back(d);
}
}
std::fclose(file);
- if (istart % 10 == 0 && istart > 0) {
+ if (istart % 1000 == 0 && istart > 0) {
std::cout << "\n";
}
std::cout << " | done " << istart << " | ";
std::vector< std::vector< RVec > > trajectory;
std::vector< RVec > reference;
+ std::vector< std::vector < std::vector < unsigned int > > > sels;
+
std::vector< int > index;
std::vector< std::string > resNms;
int frames = 0;
int tau = 0; // selectable
+
+ int tau_jump = 10;
+
int window = 0;
float crl_border = 0; // selectable
float eff_rad = 1.5; // selectable
void
SpaceTimeCorr::initAnalysis(const TrajectoryAnalysisSettings &settings, const TopologyInformation &top)
{
- ArrayRef< const int > atomind = sel_[0].atomIndices();
+ ArrayRef< const int > atomind = sel_.front().atomIndices();
index.resize(0);
for (ArrayRef< const int >::iterator ai = atomind.begin(); (ai < atomind.end()); ai++) {
index.push_back(*ai);
}
trajectory.resize(0);
+ std::string str(sel_.back().name());
+ unsigned int tempSize = std::stoul(str.substr(13, 6)) / (window / tau_jump);
+
+ sels.resize(tempSize);
+
+ for (unsigned int i = 1; i < sel_.size(); i++) {
+ std::string str(sel_[i].name());
+ atomind = sel_[i].atomIndices();
+ std::vector< unsigned int > a;
+ a.resize(0);
+ sels[std::stoul(str.substr(13, 6)) / (window / tau_jump)].push_back(a);
+ for (ArrayRef< const int >::iterator ai = atomind.begin(); (ai < atomind.end()); ai++) {
+ sels[std::stoul(str.substr(13, 6)) / (window / tau_jump)].back().push_back(*ai);
+ }
+ }
+
reference.resize(0);
if (top.hasFullTopology()) {
for (unsigned int i = 0; i < index.size(); i++) {
SpaceTimeCorr::finishAnalysis(int nframes)
{
std::vector< std::vector< std::vector< std::pair< unsigned int, unsigned int > > > > rout_new;
- unsigned int k = 1000, m = 0;
+ unsigned int k = 500, m = 0;
if (tau > 0) {
k = static_cast< unsigned int>(tau);
}
if (window == 0) {
- window = k * 2;
+ window = static_cast< int >(k * 2);
}
+ // need to double check the formula
+ //std::cout << "\nProgram requires " << 4*index.size()*index.size()*static_cast< unsigned int>(tau)*(trajectory.size() - static_cast< unsigned int>(window) - k)/tau_jump/1024/1024/1024 << " Gb of free space (may be x2)";
+
if (mode == 0) {
std::cout << "\nCorrelation's evaluation - start\n";
- correlation_evaluation(reference, trajectory, k, static_cast< unsigned int >(window), (OutPutName + "-DumpData.txt").c_str());
+ correlation_evaluation(reference, trajectory, k, static_cast< unsigned int >(window), (OutPutName + "-DumpData.txt").c_str(), tau_jump, sels);
std::cout << "Correlation's evaluation - end\n";
} else if (mode == 1) {
- rout_new.resize(trajectory.size() - static_cast< unsigned long>(window) - k);
+ rout_new.resize((trajectory.size() - static_cast< unsigned long >(window) - k)/tau_jump);
FILE *file;
file = std::fopen((OutPutName + "-DumpData.txt").c_str(), "r+");
- #pragma omp parallel for ordered schedule(dynamic) shared(rout_new) firstprivate(reference, crl_border, eff_rad, m, k)
- for(unsigned long i = 0; i < trajectory.size() - static_cast< unsigned int>(window) - k; i++) {
+ //#pragma omp parallel for ordered schedule(dynamic) shared(rout_new) firstprivate(reference, crl_border, eff_rad, m, k)
+ for(unsigned long i = 0; i < trajectory.size() - static_cast< unsigned int>(window) - k; i+= static_cast< unsigned int >(window) / tau_jump) {
std::vector< std::vector< std::vector< double > > > crltns;
std::vector< std::vector< std::pair< double, long int > > > graph;
std::vector< std::vector< unsigned int > > sub_graph;
std::vector< std::vector< std::pair< unsigned int, unsigned int > > > sub_graph_rbr;
- #pragma omp ordered
- {
+ //#pragma omp ordered
+ //{
crltns.resize(k + 1);
for (unsigned int i1 = 0; i1 < crltns.size(); i1++) {
crltns[i1].resize(trajectory.front().size());
std::cout << "\n";
}
std::cout << " | mtrxs " << i << " red | ";
- }
+ //}
graph_calculation(graph, sub_graph, sub_graph_rbr, reference, crltns, static_cast< double >(crl_border), static_cast< double >(eff_rad), m, k);
- graph_back_bone_evaluation(rout_new[i], reference.size()/*index.size()*/, graph, sub_graph, sub_graph_rbr);
+ graph_back_bone_evaluation(rout_new[i / (static_cast< unsigned int >(window) / tau_jump)], reference.size(), graph, sub_graph, sub_graph_rbr);
}
- #pragma omp barrier
+ //#pragma omp barrier
std::cout << "\n Almost - end\n";
- make_rout_file(static_cast< double >(crl_border), index, resNms, rout_new, (OutPutName + "-routs.txt").c_str(), static_cast< unsigned int >(tau), static_cast< unsigned int >(window));
+ make_rout_file(static_cast< double >(crl_border), index, resNms, rout_new, (OutPutName + "-routs.txt").c_str(), static_cast< unsigned int >(tau), static_cast< unsigned int >(window), tau_jump);
}
std::cout << "Finish Analysis - end\n\n";
}
biod.pnpi.spb.ru / alexxy / gromacs-spacetimecorr.git / commitdiff