src/spacetimecorr.cpp

   1 /*
   2  * This file is part of the GROMACS molecular simulation package.
   3  *
   4  * Copyright (c) 2013,2014,2015, by the GROMACS development team, led by
   5  * Mark Abraham, David van der Spoel, Berk Hess, and Erik Lindahl,
   6  * and including many others, as listed in the AUTHORS file in the
   7  * top-level source directory and at http://www.gromacs.org.
   8  *
   9  * GROMACS is free software; you can redistribute it and/or
  10  * modify it under the terms of the GNU Lesser General Public License
  11  * as published by the Free Software Foundation; either version 2.1
  12  * of the License, or (at your option) any later version.
  13  *
  14  * GROMACS is distributed in the hope that it will be useful,
  15  * but WITHOUT ANY WARRANTY; without even the implied warranty of
  16  * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
  17  * Lesser General Public License for more details.
  18  *
  19  * You should have received a copy of the GNU Lesser General Public
  20  * License along with GROMACS; if not, see
  21  * http://www.gnu.org/licenses, or write to the Free Software Foundation,
  22  * Inc., 51 Franklin Street, Fifth Floor, Boston, MA  02110-1301  USA.
  23  *
  24  * If you want to redistribute modifications to GROMACS, please
  25  * consider that scientific software is very special. Version
  26  * control is crucial - bugs must be traceable. We will be happy to
  27  * consider code for inclusion in the official distribution, but
  28  * derived work must not be called official GROMACS. Details are found
  29  * in the README & COPYING files - if they are missing, get the
  30  * official version at http://www.gromacs.org.
  31  *
  32  * To help us fund GROMACS development, we humbly ask that you cite
  33  * the research papers on the package. Check out http://www.gromacs.org.
  34  */
  35 /*! \internal \file
  36  * \brief
  37  * Implements gmx::analysismodules::Freevolume.
  38  *
  39  * \author Titov Anatoly <Wapuk-cobaka@yandex.ru>
  40  * \ingroup module_trajectoryanalysis
  41  */
  42
  43
  44 #include <iostream>
  45 #include <algorithm>
  46 #include <vector>
  47 #include <math.h>
  48 #include <omp.h>
  49 #include <string>
  50
  51 #include <gromacs/trajectoryanalysis.h>
  52 #include <gromacs/utility/smalloc.h>
  53 #include <gromacs/math/do_fit.h>
  54 using namespace gmx;
  55
  56 using gmx::RVec;
  57
  58 void make_correlation_matrix_file(std::vector< std::vector< std::vector< float > > > correlations, const char* file_name, int start)
  59 {
  60     FILE *file;
  61     file = std::fopen(file_name, "w+");
  62     for (int i = start; i < correlations.size(); i++) {
  63         if (correlations.size() - start > 1) {
  64             std::fprintf(file, "correlation between 'n' and 'n + %d' frames\n", i);
  65         }
  66         if (i % 50 == 0) {
  67             std::cout << "correlation between 'n' and 'n + " << i << "' frames\n";
  68         }
  69         for (int j = 0; j < correlations[i].size(); j++) {
  70             for (int f = 0; f < correlations[i][j].size(); f++) {
  71                 std::fprintf(file, "%3.4f ", correlations[i][j][f]);
  72             }
  73             std::fprintf(file, "\n");
  74         }
  75     }
  76     std::fclose(file);
  77 }
  78
  79 void make_correlation_pairs_file(std::vector< std::vector< std::vector< float > > > correlations, const char* file_name, int start)
  80 {
  81     FILE *file;
  82     file = std::fopen(file_name, "w+");
  83     for (int i = 0; i < correlations.front().size(); i++) {
  84         for (int j = 0; j < correlations.front().front().size(); j++) {
  85             std::fprintf(file, "correlation between point '%d' and point '%d'\n", i, j);
  86             for (int k = 0; k < correlations.size(); k++) {
  87                 std::fprintf(file, "%3.4f ", correlations[k][i][j]);
  88             }
  89             std::fprintf(file, "\n");
  90         }
  91         if (i % 10 == 0) {
  92             std::cout << "correlations in row " << i << " out of " << correlations.front().size() << " completed\n";
  93         }
  94     }
  95     std::fclose(file);
  96 }
  97
  98 void make_rout_file(float crl_border, std::vector< int > indx, std::vector< std::vector< std::pair< int, int > > > rout, const char* file_name)
  99 {
 100     FILE *file;
 101     file = std::fopen(file_name, "w+");
 102     std::fprintf(file, "correlations >= %0.2f\n\n", crl_border);
 103     for (int i = 0; i < rout.size(); i++) {
 104         for (int j = 0; j < rout[i].size(); j++) {
 105             std::fprintf(file, "cgo_arrow (id %3d), (id %3d), radius=0.15\n", indx[rout[i][j].first]/* + 1*/, indx[rout[i][j].second]/* + 1*/);
 106         }
 107         std::fprintf(file, "\n\n");
 108     }
 109     std::fclose(file);
 110 }
 111
 112 void make_best_corrs_graphics(std::vector< std::vector< std::vector< float > > > correlations,
 113                               std::vector< std::vector< std::pair< int, int > > > rout_pairs,
 114                               std::vector< int > indx,
 115                               const char* file_name)
 116 {
 117     FILE *file;
 118     file = std::fopen(file_name, "w+");
 119     for (int i = 0; i < rout_pairs.size(); i++) {
 120         for (int j = 0; j < rout_pairs[i].size(); j++) {
 121             std::fprintf(file, "%3d %3d\n", indx[rout_pairs[i][j].first]/* + 1*/, indx[rout_pairs[i][j].second]/* + 1*/);
 122             for (int k = 0; k < correlations.size(); k++) {
 123                 std::fprintf(file, "%3.5f ", correlations[k][rout_pairs[i][j].first][rout_pairs[i][j].second]);
 124             }
 125             std::fprintf(file, "\n");
 126         }
 127     }
 128     std::fclose(file);
 129 }
 130
 131 bool mysortfunc (std::vector< int > a, std::vector< int > b) {
 132     return (a.size() > b.size());
 133 }
 134
 135 bool isitsubset (std::vector< int > a, std::vector< int > b) {
 136     if (b.size() == 0) {
 137         return true;
 138     }
 139     std::sort(a.begin(), a.end());
 140     std::sort(b.begin(), b.end());
 141     int k = 0;
 142     for (int i = 0; i < a.size(); i++) {
 143         if (b[k] == a[i]) {
 144             k++;
 145         }
 146     }
 147     if (k == b.size()) {
 148         return true;
 149     } else {
 150         return false;
 151     }
 152 }
 153
 154 void correlation_evaluation(std::vector< std::vector< RVec > > traj, int b_frame, std::vector< std::vector< std::vector< float > > > &crl, int tauS, int tauE) {
 155     crl.resize(tauE + 1);
 156     for (int i = 0; i < crl.size(); i++) {
 157         crl[i].resize(traj.front().size());
 158         for (int j = 0; j < crl[i].size(); j++) {
 159             crl[i][j].resize(traj.front().size(), 0);
 160         }
 161     }
 162     RVec temp_zero;
 163     temp_zero[0] = 0;
 164     temp_zero[1] = 0;
 165     temp_zero[2] = 0;
 166
 167     std::vector< float > d;
 168     d.resize(traj.front().size(), 0);
 169
 170     //#pragma omp parallel shared(crl, temp_zero, d)
 171     //{
 172         //#pragma omp for schedule(dynamic)
 173         for (int i = tauS; i <= tauE; i += 1) {
 174             std::vector< RVec > medx, medy;
 175             RVec temp1, temp2;
 176             float tmp;
 177             medx.resize(traj.front().size(), temp_zero);
 178             medy.resize(traj.front().size(), temp_zero);
 179             for (int j = 0; j < traj.size() - i - 1; j++) {
 180                 for (int f = 0; f < traj.front().size(); f++) {
 181                     medx[f] += (traj[b_frame][f] - traj[j][f]);
 182                     medy[f] += (traj[b_frame][f] - traj[j + i][f]);
 183                 }
 184             }
 185             for (int j = 0; j < traj.front().size(); j++) {
 186                 tmp = traj.size() - 1;
 187                 tmp -= i;
 188                 tmp = 1 / tmp;
 189                 //tmp = 1 / (traj.size() - 1 - i);
 190                 medx[j] *= tmp;
 191                 medy[j] *= tmp;
 192             }
 193             std::vector< std::vector< float > > a, b, c;
 194             a.resize(traj.front().size(), d);
 195             b.resize(traj.front().size(), d);
 196             c.resize(traj.front().size(), d);
 197             for (int j = 0; j < traj.size() - i - 1; j++) {
 198                 for (int f1 = 0; f1 < traj.front().size(); f1++) {
 199                     for (int f2 = 0; f2 < traj.front().size(); f2++) {
 200                         temp1 = traj[b_frame][f1] - traj[j][f1] - medx[f1];
 201                         temp2 = traj[b_frame][f2] - traj[j + i][f2] - medy[f2];
 202
 203                         a[f1][f2] += (temp1[0] * temp2[0] + temp1[1] * temp2[1] + temp1[2] * temp2[2]);
 204                         b[f1][f2] += (temp1[0] * temp1[0] + temp1[1] * temp1[1] + temp1[2] * temp1[2]);
 205                         c[f1][f2] += (temp2[0] * temp2[0] + temp2[1] * temp2[1] + temp2[2] * temp2[2]);
 206                     }
 207                 }
 208             }
 209             for (int j = 0; j < traj.front().size(); j++) {
 210                 for (int f = 0; f < traj.front().size(); f++) {
 211                     crl[i][j][f] = a[j][f] / (std::sqrt(b[j][f] * c[j][f]));
 212                 }
 213             }
 214             medx.resize(0);
 215             medy.resize(0);
 216             std::cout << i << " corr done\n";
 217         }
 218     //}
 219     //#pragma omp barrier
 220 }
 221
 222 void graph_calculation(std::vector< std::vector< std::pair< float, int > > > &graph, std::vector< std::vector< int > > &s_graph, std::vector< std::vector< std::pair< int, int > > > &s_graph_rbr,
 223                       std::vector< std::vector< RVec > > traj, int b_frame,
 224                       std::vector< std::vector< std::vector< float > > > crl, float crl_b, float e_rad, int tauE) {
 225     graph.resize(traj.front().size());
 226     for (int i = 0; i < traj.front().size(); i++) {
 227         graph[i].resize(traj.front().size(), std::make_pair(0, -1));
 228     }
 229     RVec temp;
 230     for (int i = 1; i <= tauE; i++) {
 231         for (int j = 0; j < traj.front().size(); j++) {
 232             for (int f = j; f < traj.front().size(); f++) {
 233                 temp = traj[b_frame][j] - traj[b_frame][f];
 234                 if (std::max(std::abs(crl[i][j][f]), std::abs(crl[i][f][j])) >= crl_b && norm(temp) <= e_rad && std::abs(graph[j][f].first) < std::max(std::abs(crl[i][j][f]), std::abs(crl[i][f][j]))) {
 235                     if (std::abs(crl[i][j][f]) > std::abs(crl[i][f][j])) {
 236                         graph[j][f].first = crl[i][j][f];
 237                     } else {
 238                         graph[j][f].first = crl[i][f][j];
 239                     }
 240                     graph[j][f].second = i;
 241                 }
 242             }
 243         }
 244     }
 245     std::cout << "crl analysed\n";
 246     std::vector< bool > graph_flags;
 247     graph_flags.resize(traj.front().size(), true);
 248     std::vector< int > a;
 249     a.resize(0);
 250     std::vector< std::pair< int, int > > b;
 251     b.resize(0);
 252     std::vector< int > que1, que2, que3;
 253     for (int i = 0; i < traj.front().size(); i++) {
 254         if (graph_flags[i]) {
 255             s_graph.push_back(a);
 256             s_graph_rbr.push_back(b);
 257             que1.resize(0);
 258             que2.resize(0);
 259             que3.resize(0);
 260             que1.push_back(i);
 261             que3.push_back(i);
 262             graph_flags[i] = false;
 263             while(que1.size() > 0) {
 264                 que2.clear();
 265                 for (int k = 0; k < que1.size(); k++) {
 266                     for (int j = 0; j < traj.front().size(); j++) {
 267                         if (graph[que1[k]][j].second > -1 && graph_flags[j]) {
 268                             que2.push_back(j);
 269                             graph_flags[j] = false;
 270                         }
 271                     }
 272                 }
 273                 que1 = que2;
 274                 for (int j = 0; j < que2.size(); j++) {
 275                     que3.push_back(que2[j]);
 276                 }
 277             }
 278             s_graph.back() = que3;
 279             for (int j = 0; j < que3.size(); j++) {
 280                 for (int k = 0; k < traj.front().size(); k++) {
 281                     if (graph[que3[j]][k].second > -1) {
 282                         s_graph_rbr.back().push_back(std::make_pair(que3[j], k));
 283                     }
 284                 }
 285             }
 286             //std::cout << s_graph.back().size() << "   ";
 287         }
 288     }
 289 }
 290
 291 bool myfunction (const std::pair< int, float > i, const std::pair< int, float > j) {
 292     return i.second < j.second;
 293 }
 294
 295 void graph_back_bone_evaluation(std::vector< std::vector < std::pair< int, int > > > &rout_n, int indxSize,
 296                                 std::vector< std::vector< std::pair< float, int > > > graph, std::vector< std::vector< int > > s_graph, std::vector< std::vector< std::pair< int, int > > > s_graph_rbr) {
 297     std::vector< float > key;
 298     std::vector< int > path;
 299     std::vector< std::pair< int, float > > que;
 300     std::vector< std::pair< int, int > > a;
 301     int v;
 302     for (int i = 0; i < s_graph.size(); i++) {
 303         key.resize(0);
 304         path.resize(0);
 305         que.resize(0);
 306         v = 0;
 307         if (s_graph[i].size() > 2) {
 308             key.resize(indxSize, 2);
 309             path.resize(indxSize, -1);
 310             key[s_graph[i][0]] = 0;
 311             for (int j = 0; j < s_graph[i].size(); j++) {
 312                 que.push_back(std::make_pair(s_graph[i][j], key[s_graph[i][j]]));
 313             }
 314             std::sort(que.begin(), que.end(), myfunction);
 315             while (!que.empty()) {
 316                 v = que[0].first;
 317                 que.erase(que.begin());
 318                 for (int j = 0; j < s_graph_rbr[i].size(); j++) {
 319                     int u = -1;
 320                     if (s_graph_rbr[i][j].first == v) {
 321                         u = s_graph_rbr[i][j].second;
 322                     } else if (s_graph_rbr[i][j].second == v) {
 323                         u = s_graph_rbr[i][j].first;
 324                     }
 325                     bool flag = false;
 326                     int pos = -1;
 327                     for (int k = 0; k < que.size(); k++) {
 328                         if (que[k].first == u) {
 329                             flag = true;
 330                             pos = k;
 331                             k = que.size() + 1;
 332                         }
 333                     }
 334                     if (flag && key[u] > 1 - std::abs(graph[v][u].first)) {
 335                         path[u] = v;
 336                         key[u] = 1 - std::abs(graph[v][u].first);
 337                         que[pos].second = key[u];
 338                         sort(que.begin(), que.end(), myfunction);
 339                     }
 340                 }
 341             }
 342             a.resize(0);
 343             rout_n.push_back(a);
 344             for (int j = 0; j < indxSize; j++) {
 345                 if (path[j] != -1) {
 346                     rout_n.back().push_back(std::make_pair(j, path[j]));
 347                 }
 348             }
 349         }
 350     }
 351 }
 352
 353 /*! \brief
 354  * \ingroup module_trajectoryanalysis
 355  */
 356
 357 class SpaceTimeCorr : public TrajectoryAnalysisModule
 358 {
 359     public:
 360
 361         SpaceTimeCorr();
 362         virtual ~SpaceTimeCorr();
 363
 364         //! Set the options and setting
 365         virtual void initOptions(IOptionsContainer          *options,
 366                                  TrajectoryAnalysisSettings *settings);
 367
 368         //! First routine called by the analysis framework
 369         // virtual void initAnalysis(const t_trxframe &fr, t_pbc *pbc);
 370         virtual void initAnalysis(const TrajectoryAnalysisSettings &settings,
 371                                   const TopologyInformation        &top);
 372
 373         virtual void initAfterFirstFrame(const TrajectoryAnalysisSettings   &settings,
 374                                          const t_trxframe                   &fr);
 375
 376         //! Call for each frame of the trajectory
 377         // virtual void analyzeFrame(const t_trxframe &fr, t_pbc *pbc);
 378         virtual void analyzeFrame(int frnr, const t_trxframe &fr, t_pbc *pbc,
 379                                   TrajectoryAnalysisModuleData *pdata);
 380
 381         //! Last routine called by the analysis framework
 382         // virtual void finishAnalysis(t_pbc *pbc);
 383         virtual void finishAnalysis(int nframes);
 384
 385         //! Routine to write output, that is additional over the built-in
 386         virtual void writeOutput();
 387
 388     private:
 389
 390         SelectionList                                               sel_;
 391
 392         std::vector< std::vector< RVec > >                          trajectory;
 393
 394         std::vector< int >                                          index;
 395         int                                                         frames              = 0;
 396         int                                                         basic_frame         = 0;
 397         int                                                         tau                 = 0; // selectable
 398         float                                                       crl_border          = 0; // selectable
 399         float                                                       eff_rad             = 1.5; // selectable
 400         std::string                                                 OutPutName; // selectable
 401         // Copy and assign disallowed by base.
 402 };
 403
 404 SpaceTimeCorr::SpaceTimeCorr(): TrajectoryAnalysisModule()
 405 {
 406 }
 407
 408 SpaceTimeCorr::~SpaceTimeCorr()
 409 {
 410 }
 411
 412 void
 413 SpaceTimeCorr::initOptions(IOptionsContainer          *options,
 414                      TrajectoryAnalysisSettings *settings)
 415 {
 416     static const char *const desc[] = {
 417         "[THISMODULE] to be done"
 418     };
 419     // Add the descriptive text (program help text) to the options
 420     settings->setHelpText(desc);
 421     // Add option for output file name
 422     //options->addOption(FileNameOption("on").filetype(eftIndex).outputFile()
 423     //                        .store(&fnNdx_).defaultBasename("domains")
 424     //                        .description("Index file from the domains"));
 425     // Add option for tau constant
 426     options->addOption(gmx::IntegerOption("tau")
 427                             .store(&tau)
 428                             .description("number of frames for time travel"));
 429     // Add option for crl_border constant
 430     options->addOption(FloatOption("crl")
 431                             .store(&crl_border)
 432                             .description("make graph based on corrs > constant"));
 433     // Add option for eff_rad constant
 434     options->addOption(FloatOption("ef_rad")
 435                             .store(&eff_rad)
 436                             .description("effective radius for atoms to evaluate corrs"));
 437     // Add option for selection list
 438     options->addOption(SelectionOption("select_domains_and_residue").storeVector(&sel_)
 439                            .required().dynamicMask().multiValue()
 440                            .description("Domains to form rigid skeleton"));
 441     options->addOption(StringOption("out_put")
 442                             .store(&OutPutName)
 443                             .description("<your name here> + <local file tag>.txt"));
 444     // Control input settings
 445     settings->setFlags(TrajectoryAnalysisSettings::efNoUserPBC);
 446     settings->setPBC(true);
 447 }
 448
 449 void
 450 SpaceTimeCorr::initAnalysis(const TrajectoryAnalysisSettings &settings,
 451                       const TopologyInformation        &top)
 452 {
 453     ArrayRef< const int > atomind = sel_[0].atomIndices();
 454     index.resize(0);
 455     for (ArrayRef< const int >::iterator ai = atomind.begin(); (ai < atomind.end()); ai++) {
 456         index.push_back(*ai);
 457     }
 458     trajectory.resize(0);
 459 }
 460
 461 void
 462 SpaceTimeCorr::initAfterFirstFrame(const TrajectoryAnalysisSettings       &settings,
 463                              const t_trxframe                       &fr)
 464 {
 465 }
 466
 467 // -s '/home/toluk/Develop/samples/reca_rd/reca_rd.mono.xtc' -f '/home/toluk/Develop/samples/reca_rd/reca_rd.mono.xtc' -n '/home/toluk/Develop/samples/reca_rd/test6.ndx' -sf '/home/toluk/Develop/samples/reca_rd/SelectionList5'  -tau 5 -crl 0.10
 468 // -s '/home/toluk/Develop/samples/reca_rd/reca_rd.mono.tpr' -f '/home/toluk/Develop/samples/reca_rd/reca_rd.mono.xtc' -n '/home/toluk/Develop/samples/reca_rd/CorrsTestDomainsNfit.ndx' -sf '/home/toluk/Develop/samples/reca_rd/SelectionListDomainsNFit' -tau 5000 -crl 0.75 -ef_rad 1
 469 // -s '/home/toluk/Develop/samples/reca_rd/reca_rd.mono.tpr' -f '/home/toluk/Develop/samples/reca_rd/reca_rd.mono.xtc' -n '/home/toluk/Develop/samples/reca_rd/TestCa.ndx' -sf '/home/toluk/Develop/samples/reca_rd/SelListCa' -tau 100 -crl 0.75 -ef_rad 1
 470 // -s '/home/toluk/Develop/samples/reca_rd/reca_rd.mono.tpr' -f '/home/toluk/Develop/samples/reca_rd/reca_rd.mono.xtc' -n '/home/toluk/Develop/samples/reca_rd/testCa.ndx' -sf '/home/toluk/Develop/samples/reca_rd/SelListCa' -tau 9000 -crl 0.60 -ef_rad 1
 471
 472 // -s '/home/toluk/Develop/samples/CubeTermal/CUBETermalTest1.pdb' -f '/home/toluk/Develop/samples/CubeTermal/CUBETermalTest.pdb' -n '/home/toluk/Develop/samples/CubeTermal/testCa.ndx' -sf '/home/toluk/Develop/samples/CubeTermal/SelListCa' -tau 900 -crl 0.20 -ef_rad 1
 473
 474 void
 475 SpaceTimeCorr::analyzeFrame(int frnr, const t_trxframe &fr, t_pbc *pbc,
 476                       TrajectoryAnalysisModuleData *pdata)
 477 {
 478     trajectory.resize(trajectory.size() + 1);
 479     trajectory.back().resize(index.size());
 480     for (int i = 0; i < index.size(); i++) {
 481         trajectory.back()[i] = fr.x[index[i]];
 482     }
 483     frames++;
 484 }
 485
 486 void
 487 SpaceTimeCorr::finishAnalysis(int nframes)
 488 {
 489     std::vector< std::vector< std::vector< float > > > crltns;
 490     std::vector< std::vector< std::pair< float, int > > > graph;
 491     std::vector< std::vector< int > > sub_graph;
 492     std::vector< std::vector< std::pair< int, int > > > sub_graph_rbr, rout_new;
 493     int k = 1000, m = 0;
 494     if (tau > 0) {
 495         k = tau;
 496     }
 497
 498     std::cout << "\nCorrelation's evaluation - start\n";
 499
 500     correlation_evaluation(trajectory, basic_frame, crltns, m, k);
 501
 502     std::cout << "Correlation's evaluation - end\n" << "graph evaluation: start\n";
 503
 504     graph_calculation(graph, sub_graph, sub_graph_rbr, trajectory, basic_frame, crltns, crl_border, eff_rad, k);
 505
 506     std::cout << "graph evaluation: end\n" << "routs evaluation: start\n";
 507
 508     graph_back_bone_evaluation(rout_new, index.size(), graph, sub_graph, sub_graph_rbr);
 509
 510     std::cout << "routs evaluation: end\n";
 511
 512     make_correlation_matrix_file(crltns, (OutPutName + "_matrix.txt").c_str(), 0);
 513     std::cout << "corelation matrix printed\n";
 514     make_correlation_pairs_file(crltns, (OutPutName + "_pairs.txt").c_str(), 0);
 515     std::cout << "corelation pairs printed\n";
 516     make_rout_file(crl_border, index, rout_new, (OutPutName + "_routs.txt").c_str());
 517     std::cout << "corelation routs printed\n";
 518     make_best_corrs_graphics(crltns, rout_new, index, (OutPutName + "_routs_graphics.txt").c_str());
 519     std::cout << "corelation routs' pairs' graphics printed\n";
 520
 521     std::cout << "Finish Analysis - end\n\n";
 522 }
 523
 524 void
 525 SpaceTimeCorr::writeOutput()
 526 {
 527
 528 }
 529
 530 /*! \brief
 531  * The main function for the analysis template.
 532  */
 533 int
 534 main(int argc, char *argv[])
 535 {
 536     return gmx::TrajectoryAnalysisCommandLineRunner::runAsMain<SpaceTimeCorr>(argc, argv);
 537 }