From 5b4d376fe363d7499f2f9d12bb58f9b9726a181a Mon Sep 17 00:00:00 2001 From: Anatoly Date: Wed, 20 Nov 2019 16:59:40 +0300 Subject: [PATCH] - added trj fitting --- src/spacetimecorr.cpp | 265 +++++++++++++++++++++++++++++++++++++----- 1 file changed, 239 insertions(+), 26 deletions(-) diff --git a/src/spacetimecorr.cpp b/src/spacetimecorr.cpp index a0ad64f..8eed275 100644 --- a/src/spacetimecorr.cpp +++ b/src/spacetimecorr.cpp @@ -49,6 +49,7 @@ #include #include #include +#include #include #include @@ -61,6 +62,121 @@ using namespace gmx; 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; @@ -71,7 +187,7 @@ const std::vector< std::vector< std::vector< double > > > read_correlation_matri 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); @@ -91,7 +207,7 @@ const std::vector< std::vector< std::vector< double > > > read_correlation_matri } 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; @@ -103,7 +219,7 @@ unsigned int posSrch(std::vector< T > a, T b) { 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+"); @@ -121,9 +237,9 @@ void make_rout_file(double crl_border, std::vector< int > indx, std::vector< std } } 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++) { @@ -134,23 +250,27 @@ void make_rout_file(double crl_border, std::vector< int > indx, std::vector< std 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)); @@ -222,11 +342,24 @@ bool isitsubset (std::vector< int > a, std::vector< int > b) { } } -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++) { @@ -234,6 +367,59 @@ void correlation_evaluation(const std::vector< RVec > ref, const std::vector< st 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()); @@ -249,6 +435,9 @@ void correlation_evaluation(const std::vector< RVec > ref, const std::vector< st 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); @@ -305,7 +494,7 @@ void correlation_evaluation(const std::vector< RVec > ref, const std::vector< st } } std::fclose(file); - if (istart % 10 == 0 && istart > 0) { + if (istart % 1000 == 0 && istart > 0) { std::cout << "\n"; } std::cout << " | done " << istart << " | "; @@ -502,11 +691,16 @@ class SpaceTimeCorr : public TrajectoryAnalysisModule 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 @@ -573,7 +767,7 @@ SpaceTimeCorr::initOptions(IOptionsContainer *options, 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); @@ -581,6 +775,22 @@ SpaceTimeCorr::initAnalysis(const TrajectoryAnalysisSettings &settings, const To } 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++) { @@ -611,33 +821,36 @@ void 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()); @@ -660,15 +873,15 @@ SpaceTimeCorr::finishAnalysis(int nframes) 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"; } -- 2.22.0